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THE STAR FORMATION RATE OF TURBULENT MAGNETIZED CLOUDS: 
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ABSTRACT 

The role of turbulence and magnetic fields is studied for star formation in molecular clouds. We 
derive and compare six theoretical models for the star formation rate (SFR) — the Krumholz & McKee 
(KM), Padoan & Nordlund (PN), and Hennebelle & Chabrier (HC) models, and three multi-freefall 
versions of these, suggested by HC — all based on integrals over the log-normal distribution of turbulent 
gas. We extend all theories to include magnetic fields, and show that the SFR depends on four basic 
parameters: (1) virial parameter a v ir! (2) sonic Mach number Ai\ (3) turbulent forcing parameter 
b, which is a measure for the fraction of energy driven in compressive modes; and (4) plasma /3 — 
2Ai\/A4 2 with the Alfven Mach number A4a- We compare all six theories with MHD simulations, 
covering cloud masses of 300 to 4 x 10 6 M and Mach numbers M. = 3-50 and A4a = l-oo, with 
solenoidal (6 = 1/3), mixed (b — 0.4) and compressive turbulent (b — 1) forcings. We find that the 
SFR increases by a factor of four between M = 5 and 50 for compressive turbulent forcing and a V ir ~ 1- 
Comparing forcing parameters, we see that the SFR is more than 10 x higher with compressive than 
solenoidal forcing for A4 = 10 simulations. The SFR and fragmentation are both reduced by a factor 
of two in strongly magnetized, trans- Alfvenic turbulence compared to hydrodynamic turbulence. All 
simulations are fit simultaneously by the multi-freefall KM and multi-freefall PN theories within a 
factor of two over two orders of magnitude in SFR. The simulated SFRs cover the range and correlation 
of SFR column density with gas column density observed in Galactic clouds, and agree well for star 
formation efficiencies SFE = 1%— 10% and local efficiencies e = 0.3-0.7 due to feedback. We conclude 
that the SFR is primarily controlled by interstellar turbulence, with a secondary effect coming from 
magnetic fields. 

Subject headings: ISM: clouds - ISM: kinematics and dynamics - ISM: structure - Magnetohydrody- 
namics (MHD) - Stars: formation - Turbulence 
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1. INTRODUCTION 

Stars form in turbulent, magnetized molecular clouds, 
as observed in the Milky Way and in other galaxies. Yet, 
the basic physical processes controlling star formation 
are still poorly understood. Observations of star-forming 
clouds show that the star formation rate (SFR) column 
density Ssfr varies over four orders of magnitude and 
exhibits a positive correlation w ith the gas surface den- 
sity Sgas ( Heiderman et al.|2010 ), suggesting that denser 
gas forms stars at a higher rate. This engenders the 
central question of how the gas is locally compressed 
in the interstellar medium, such that dense cores can 
form and eventually become unstable under their own 
gravitational attraction to form stars. Gas compression 
in shocks, induced by large-scale supersonic turbulence 
might be a key — if not the key process — setting the ini- 
tia l conditions for star formation (see, e.g., the review s 
b y [Mac Low fc Klessen||2004} |McKee fc Ostriker|[2007l ) . 

Based on molecular cloua masses m the range 
M c = 100 to 10 7 M Q and temperatures T < 20 K, the 
clouds should be highly Jeans-unstable and would thus 
collapse globally. However, molecular clouds do not 
show systematic, global collapse motions. If they 
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did, the average Galac tic SFR in the Milky Way. 

l^M^yr" 1 (|Robitaille fc Whitneyj |2010 



SFR 



MW 



|Chomiuk & Povich 2011) would be about t wo~oTders o: 
magn itude high e r than the observed yalue ( Zuckerman 
& Palmerl 1974 Zuckerman & Evans 1974)! However, 



this stability analysis only takes thermal pressure into 
account. In reality, clouds ar e magnetized and subject 



to strong turbulent moti ons (Scalo & Elmegreen 2004 
Elmegreen fc ScalopOOlt 

Originally, it has been thought that primarily magnetic 
fields would provide stability against fast global collapse, 
and that only after the neutral species have slowly dif- 
fused through the charged particles, star formatio n would 
occur in the central regions of magnetized clo uds ( |Mestel| 
|fc Spitzer|1956} |Mouschovias|1976[ [Shu|1983[ ) . In this so- 
called ambipolar-diiiusion process, magnetic flux is left 
behind in the envelope, while the mass increases in the 
cloud core. Thus, star formation regulated by ambipo- 
lar diffusion predicts a higher mass-to-flux ratio in the 
cores than in the envelopes of th e clouds, which is - 
however — typically not observed ( Crutcher et al."||2009 
Mouschovias fc Tassis|2009| |Lunttila et al |2009| |Santos^ 
Lima et al.||2010| |Lazarian et al. 2012 Bertram et al. 



2012) 



~Kn alternative scenario is that the observed supersonic 
random motio ns flZuckerman fc Palmer|1974[|Zuckerman 
fc Evans|1974[|Larson|19caHSolomon et, al.|1987j |Falgar' 
one et a'.||1992| |Ossenkopf fc Mac Low||2002| |Heyer fc 



Brunt 12"0~04| |Schneider et al. 201 1[ |Roman-Duval et al 
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20111 regulate star formation. In this picture, turbulent 



energy stabilizes the clouds on large scales, but at the 
same time, supersonic turbulence induces local compres- 
sions, producing filaments and cores, which are the pro- 
genitors of stars. Eventually, both turbulence and mag- 
netic fields play their parts; the only question is: which 
one is the dominant controlling factor of star formation? 

The aim of this paper is to advance our understanding 
of the relevant physical processes and their parameters 
controlling the conversion of dense gas into stars, and to 
explain the observed variations of the SFR column den- 
sity. We develop and compare six predictive theories — 
the original Krumholz & McKee (KM), Padoan & Nord- 
lund (PN), and Hennebelle & Chabrier (HC) theories, 
and multi-freefall versions of theses three — , which are 
all based on integrals over the turbulent density proba- 
bility distribution function (PDF), explained in detail in 
the next section. We extend the KM and HC theories, as 
well as all the multi-freefall theories to include magnetic 
fields. We evaluate the relative importance of turbulence, 
its forcing characteristics, and magnetic fields in control- 
ling the SFR and show that the SFR depends on the 
following four basic parameters: 



1. virial parameter a V j r = 2E\ c \ n /\E l 

2. sonic Mach number A4 = <7v/c s 



grav 1 1 



3. turbulent forcing parameter b, with purely 
solenoidal (divergence-free) forcing parameterized 
by b — 1/3, mixed forcing by b — 0.4 and purely 
compressive (curl- free) forcing by b = 1, and 

4. the ratio of thermal to magnetic pressure f3 = 
2M\/M 2 with the Alfven Mach number Ma- 

We test all six theories with numerical simulations of 
supersonic, magnetized turbulence including self-gravity 
and sink particles to capture dense, collapsing, star- 
forming gas. We find that the multi-freefall KM and PN 
models including magnetic fields provide the best fits to 
our numerical simulations with typical uncertainties of 
less than a factor of two. This is an encouraging agree- 
ment, given that the SFR varies by two orders of mag- 
nitude in the simulations, depending on the four basic 
cloud parameters listed above. 

Comparing our numerical experiments with SFRs mea- 
sured in Galactic star-forming regions, we find that for 
typical star formation efficiencies of SFE = 1%-10%, the 
best-fit local efficiencies due to radiative and mechanical 
feedback from jets, winds, expanding shells or outflows 
driven by young stellar objects are e = 0.3-0.7 with a 
best-fit value of e w 0.5 for SFE = 3%. This suggests 
that a fraction e « 0.3-0.7 of all the infalling gas onto 
a typical protostellar core is accreted by the protostar, 
while a fraction (1 — e) ps 0.3-0.7 is re-injected into the 
interstellar medium by jets, winds, and outflows. We 
find good agreement between the numerical simulations 
and Galactic observations, suggesting that the observed 
variations in £sfr with E gas are a result of different com- 
binations of the four basic parameters controlling the 
SFR: a v i r , M, b, and /?, as listed above. Since molec- 
ular clouds are often characterized by virial parameters 
of order unity, we conclude that the degree of compres- 
sion induced by the turbulent forcing and sonic Mach 



number have the strongest influence on the SFR, caus- 
ing variations by more than an order of magnitude, while 
magnetic fields can account for reductions of the SFR by 
a factor of two. 

In Section [2l we introduce and discuss the six ana- 
lytic theories for the SFR, based on the turbulent density 
PDF, derive and discuss their dependencies, add mag- 
netic fields to the theories that did not include magnetic- 
field effects in previous derivations, and compare them 
with each other in detail. We then test the analytic 
theories with numerical simulations of supersonic, mag- 
netized turbulence by varying the sonic Mach number 
(A4 = 3-50), the forcing of the turbulence (solenoidal, 
mixed, compressive), and the magnetic- field strength 
(yielding Alfven Mach numbers A4a = 1.3-oo) to cover 
a comprehensive range of cloud parameters. The simu- 
lation methods and setups are explained in Section [3j 
A detailed time-evolution analysis of column density, 
magnetic-field morphology, and fragmentation proper- 
ties is presented in Section [4] In Section [5j we com- 
pare the SFRs measured in the magnetohydrodynamic 
(MHD) simulations with the six theoretical models, and 
determine the best-fit theory parameters that are uni- 
versally applicable and fit all our simulations simultane- 
ously Section [6] presents a comparison of SFR column 
densities in the simulations with observations of Galac- 
tic clouds. We discuss limitations of the theoretical and 
numerical models, as well as limitations in the compar- 
ison with observations in Section [7J Finally, we list our 
conclusions and summarize the most important results 
in Section |8l Here, we study the SF R in detail, while in 
Paper II ( |Federrath fc Klessen 2012), we concentrate on 
the star formation efficiency (SFE). 

2. THE SFR FROM THE STATISTICS OF SUPERSONIC 
MAGNETIZED TURBULENCE 

2.1. The Density PDF 

The probability density function (PDF) of the gas den- 
sity in a turbulent medium — such as a molecular cloud — 
is the key ingredient for analytic models of star forma- 
tion. A log-normal density PDF has been used to ex- 
plain the mass distribution of cores and stars, the core 
mass f unction (CMF) and the stellar initial mass function 
(IMF ) (|Padoan fc Nordlund|2002[ [Hennebelle fc Chabrier 
200 9i|Elmegreen|2011|IVeltchev et al.|201H|Donkov 



2008, 



et al 



2012 Parravano et al. 2012: Hopkins 2012), the 



Kennlcutt-Schmidt rela tion (Kr umholz fc JVLcKeef 2005 
l Tass^|20Q7|, the SFE ([ Ehncgrcen 2008) , and tTTe^FTi 
flKrumholz fc McKee|| 2 p 05| JPadoan fc Nordlundl^Oll, 
|Hennebelle fc Chabrier 201 Here we concentrate on 
the SFR and derive its basic dependencies. 
The log-normal PDF of the gas density is defined as, 



1 ( (s-s ) 2 

p s (s) = / _ — - exp 



expressed in terms of the logarithmic density 
s = In (p/po) . 



(1) 



(2) 



The PDF is a normal (Gaussian) distribution in s, mean- 
ing it is a log-normal distribution in p. The quantities 
Po and so denote the mean density and mean logarith- 
mic density, the latter of which is related to the standard 
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deviation a s by 



so 



(3) 



due to the normalization and mass-conservation con- 



straints of th e PDF ( Vazquez-Semadeni||1994 Federrath 
et al. 2008b). The reason to use s instead of p in the 
context ot the density PDF, is that s is dimensionlcss, 
and that the PDF of s is Gaussian unlike the PDF of p. 
This is because the distribution of p is generated by a 
multiplicative process in which shocks are amplified by 
other shocks as they collide and interact in isothermal 
supersonic turbulence, with the loca l Mach number be- 
ing independent of the local density (|Vazquez-Semadeni 
Tassot k Vazquez-Semadeni"||1998| |Kritsuk et al. 
Federrath et al.|2010b| ). Since s oc m(p) as defined 



1994 



2007 



in Equation Q , this multiplicative process in p turns into 
an additive process in s. Following the central limit the- 
orem, a large sum of random variables produces a Gaus- 
sian distribution, and thus only p s is Gaussian, while p p 
is not. However, p s can be easily transfor med into p be - 
cause p s ds = p p dp, and thus p p — p s /p ( |Li et al.||2003 ). 
We will omit the index s in p s in the following and simply 
use p(s) for the PDF given by Equation ([lj. 

As soon as significant collapse sets in, the dens ity PDF 
develops a power-law tail at h igh densities (e.g., Klessen 
|2000[ |Kainulainen et aL 2009[ ) , which is discuss ed m more 
detail i n Section 7.1.1 below, and in Paper II ( Federrath 



k Klessen||20l2"l ) 



2.2. The Standard Deviation of Density Fluctuations in 
Supersonic, Magnetized Turbulence 

The standard deviation o~ s in Equation ([I]), which is 
a measure of how much the density varies m a turbu- 
lent medium, depends on (1) the amount of compres- 
sion induced by the turbulent forcing mechanism, (2) 
the Mach number, and (3) the degree of magnetization. 
First, the turbulent energy injection mechanism deter- 
mines the amount of compression induced directly by 
driving turbulence in the interstellar medium (ISM). Var- 
ious turbulent dri ving mechanisms have been discussed 
and compared in Mac Low k Klessen| (120041) . For in- 
stance, expanding supernov a shells"! Balsara et aX1|200 4 
de Avillez k Breitschwerdt||2005| |Tamburro et al.l |200! 

sir 



of stars dMcKee 1989; Krumholz et al.|2006[|Gritschneder 


et al.||2009| |Peters et al. |2010| 


Goldbaum et al.||2011 


) as 



qp: 

(Elme green||20'09"l) and gravitational contraction flHoyle 
1953; Vazquez-Scmadeni et al.l 11998 Klessen k Hen- 



531 
nebel 



le 2010; Elm. 
:! are likely e 



|Kle 



egreen k Burkert|2010 Federrath et a 



201 lep are likely exciting a considerable amount of com 



pressible modes that will directly lead to compression, 
and thus to higher density contrasts on molecular cloud 
scales in the ISM, while g alactic rotation and magnetoro- 
tational instabilities (e.g., |Piontek k Ostriker|2004|[2007 1 
are likely producing more solenoidal modes, second, 
higher Mach numbers M. lead to stronger shocks and 
thus to higher density contrasts. For instance, the den- 
sity jump in a non-magnetized, isothermal shock is pro- 
portional to M 2 . Finally, higher magnetization dampens 
density fluctuations as magnetic fields ac t like a cushion 
due to the additional m agnetic pressure (Ostriker et al 
2001[|Price et al.||2011[ ). 



The actual dependence of turbulent density fluctua- 
tions a s on the three parameters above (forcing, Mach 
number, and magnetic field) can be derived from the 
shock jump conditions of an individual MHD shock, and 
then averaged over a who le ensemble of s uch s hocks 
( |Padoan k Nordlund|[20TT| ). |Molina et al] ( |2012[ ) pro- 
vide a rigorous derivation of a s for different correlations 
of the magnetic field with density. They distinguish three 
cases, B oc p°, B oc p 1//2 , and B oc p 1 . For the interme- 
diate case, Molina et al. (2012) derive 



= In 1 + b 2 M' 



13 



/3 + 1 



(4) 



which is s imilar to the relation derived in IPadoan 13 
Nordlund (2011), except for the factor b 2 , expla ined be- 
low, andexcept for the definition of (3, for which [Padoan| 
k Nordlund (2011 ) only take post-shock gas into account 
(see the m ore extended discussion on this issue in Sec- 
The case B oc p 1 is similar to the interme- 



2.4.21 



tion 

diate case, but is a rather extreme MHD case because 
magnetic-field lines are assumed to be oriented only per- 
pendicular to the flow direction. So is the other extreme 
case in which the magnetic field is assumed to be par- 
allel to the flow, yielding B oc p°. In the more realis- 
tic case of turbulent flows, field lines become tangled, 
and the B-p correlation is a combination of compres- 
sion of field lines and turbulent dynamo amplification 



( Schleicher et al. 12010 



20TTc ' 'lurk et al 2012 



Sur et al. 2010 Federrath et al 



[Schober et al.|26l2[ ). In a three 
;h a random distribution of flov 



dimensional system with a random distribution of flow 
velocities and magnetic-field orientations, B oc p 1 ! 2 pro- 
vides a reasonable intermediate dependence. We will 
thus only co nsider B oc p 1//2 here, which is favored by 



simulations (Padoan fc Nordlund 1999 Collins et al 



2011 |Molina~et al.||2012| , and also close to what is sug- 
gested from observations of magnetic fields in molecular 



clouds ( |Crutcher et~al]|2010| f 



In the case of B oc /9 U , i.e., for no density correla- 
tion of the magnetic field, Equation (|4l reduces to the 



pression, a 2 = In (l + b 2 M 2 ) with j3 — > oo (e.g., 


Padoan 


et al. 


1997; Passot & Vazquez-Semadeni| 1998; 


Jstriker 


et al. 


2001; Lcmaster k Stone||2008| |Federrath et al. 


20081) 


| |Price et al.l|201ip as a necessary condition in the 



tion Bl) are the turbulent forcing parameter, the rms 
sonic Mach number, and the ratio of thermal to magnetic 
pressure, plasma j3 = i"th/-Pmag- Using the definitions of 
the thermal pressure for an isothermal equation of state 



Pth = pc-l, magnetic pressure P n 



B 2 /(8n), Alfven 



velocity vj^ = B /(Airp), sonic and Alfven Mach num- 
bers, A4 = <tv/ c s and A4a — ov/^a, the plasma beta 
can be expressed as j3 = 2c 2 /v\ = 2M\/ M 2 . These 
are all dimensionless numbers, rendering them particu- 
larly useful because they determine the basic properties 
of turbulent plasmas and can thus be compared directly 
for any such system. Equation Q can thus also be writ- 



The observationally determined expo nent of the B 



lation is quite uncert ain. 



( |1999| find B oc p°- 



p corre- 
, while 



Crutcher 

|Crutcher et al.| |2010| find B oc p u below gas densities of 300 cm 
and B oc p u ot> above. For simplicity, we adopt Equation Q, de 
rived for the intermediate case, B oc p 1 / 2 . 



4 
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ten as 



a* = In 1 + b 2 M' 



2M\ 



M 2 



2M\ 



(5) 



The forcing parameter b was shown to vary smoothly 
between b ss 1/3 for purely solenoidal (divergence- free) 
forcing, and fowl for purely compressive (curly- free) 



et al. 


2009; 


Federrath et al. 


lisoiobiiy 


Micic et al. 


2012| |Konstanc 


in et al. 



Seifricd et al. 2011b; 



2012a). 



A stochas- 



leads to b ss 0.4 (see Figure 8 in i Federrath et al.|| 2 010b|) . 

Using numerical simulations, Molina et al.| Q2012| ) 
found that Equations Q and |5| work well in the regime 
M.A <; 2, while for M.a ^ 2, the assumption of isotropy 
entering the analytic derivation of Equations Q and (J5J) 
breaks down, so we only apply them in the super- Alfvemc 
regime in all the following. 

2.3. Basics of the SFR Derivation 

Here we present an analytic derivation of the SFR from 
the statistics of supersonic, isothermal, magnetized tur- 
bulence. The main ingredient for this analytic deriva- 
tion is an integral over the density PDF, Equation ([I]), 
in order to estimate the gas mass above a given density 
threshold, contributing to star formation. We will com- 
pare different ways of estimating the density threshold, 
which is the main difference between the thre e most suc- 
cessful, existing analytic models for the SFR flKrumholz 



fe McKee|2 005; Pa doan fc Nordlund|201l| |Heimebelle E 
Chabrier|201ip . We will express all quantities in terms of 
dimensioniess numbers, in order to simplify the deriva- 
tion and to make it more general. We follow the stan- 
dard terminology and use the Star Formation Rate per 
Free fa ll Time (SFRgf), as coined by Krumholz & McKee 
(2005), which is the mass fraction going into stars per 



time, where the time is expressed in units of the mean 
f reef all time. 

The SFR in units of M yr _1 can be computed by 
scaling SFRff with the real cloud mass M c and the actual 
freefall time evaluated at the mean density of the cloud, 
tfi(po)-- 

SFR ee -4^t SFR ff . (6) 
tff(po) 

Note that this definition of SFRff is different from the 



definition us ed in Krumholz & Tan (2007) and [Krumholz 
et al. (20121, who use freefall times estimated at differ- 



ent densities and/or use a definition based on column 
densities, such that the values of SFRff quoted in those 
studies and the ones computed here cannot be directly 
compared. For instance, given an SFR for fixed M c , the 
dimensioniess value of SFRff would be much smaller, if 
the freefall time at a high-density tracer was used rather 
than the freefall time at the mean density of the cloud 
because tg(p > po) is shorter than tg(po). 

The basic idea for an analytic model of SFRff is to 
integrate the log- normal density PDF, Equation (fll), 
weighted by p/po to get the mass fraction of gas with 
density above a critic al density s cr ft (to be determined 
below in Section |2.4[ ), and weighted by a freefall-time 
factor to construct a dimensioniess mass rate: 



SFRff = — 



ts(po) P 
_ ts(p) Po 



p(s) ds . 



(7) 



Note that the factor tg (po)Aff(p) appears inside the in- 
tegral because gas with different densities has different 
freefall times, 



tsip) = (sib) 



1/2 



(8) 



which sho uld be taken into account in the most general 
case (see Hennebelle & Chabrier 2011). Previous esti- 
mates for SFRff either used a factor te(po)/ts(po) = 1 
(Krumholz & McKee 2005) , or a factor tff (po)/tff (p C rit ) 



with p cr it = po exp (s cr j t ) ( Padoan fc Nordlund||2011 ) 
both of which are independent of density and were thus 
pulled out of the integral. We will show, however, that 
it is crucial to take the multi-freefall nature of gas with 
different densities into account to obtain better models 
for SFRff. 

The constant factor e in Equation ¥f\ accounts for the 
fact that only a certain fraction of the gas above s cr j t 
might actually go into stars. Since individual stars form 
in accretion disks from which powerful jets, winds, and 
outflows are launched during the process of stellar birth, 
it is likely that a certain fraction of the accreted material 
is re-injected into the ISM, thus leading to e < 1. The- 
oretical upper limits are in the range e « 0.25-0.7 (e.g., 
|Matzner fc Mc Kee 2000]). The observed d isplacem ent of 
the characteristic mass in the IMF (e.g., |Kroupa||200Tl 
Chabrier 2003) with respect to the CMF (e.g., |Johnst = one] 
et al. ' 2000 jna s been taken to argue that e~ migh t be 
aroundT J^-O .5 j Alves et al |2007[ | Andre et al.|2010[ ) ; see 
however | Ward et ai.| ( 2012 ) . 



The factor l/<t> t m Equation ^ is also of order unity 
and accounts for the uncertainty in t he timescale factor 



tff(pp)/tf{( p), originally introduced in Krumholz & Mc 
|Kee| p005 ) . We will determine the best-fit values of e and 
l/(f>t in Sections [4] and [6j when we compare the theories 
with simulations and observations. 

2.4. Six Models for the SFR 

In the following, we will solve Equation Q, using dif- 
ferent density thresholds s cr i t , accor ding to the previ- 
ous analyt ic stu dies of the SFR by |Krumholz fc Mc-| 



Keel ([20051 KM), |Padoan fc Nordlund| Q2011[ PJN),~and 
Hennebelle fc Chabrier| fl2011| HCf | We distinguish six 
cases, named l KM\ 'PN'THTT, and%iulti-ff KM', 'multi- 



ff PN', 'multi-ff HC as distinguished in |Hennebelle &1 
Chabrier (2011). The fir st three represent th e origi - 

" 20051) 



nal analytic derivations by Krumholz fc McKee| , . 

Padoan fc Nordlund (2011), and Hennebelle fc Chabrier| 



( |201l| , while the last set of three are all based on the 
rnuiti-freefall expression of the integral Q. The dif- 
ference for this last set of three is only the model for 
the critical density, i.e., the lower limit of the integral. 
We note that the ideas inherent in each of the original 



turbulence-regulated SFR. Krum holz fc McKee 



developed the basic model, .Padoan & INordlund 





2005 




2U11 



extended it to inc lude magnetic fields, and |Henne belle 



& Chabrier ( |2011[ ) improved all models by introducing 



4 Note that the critical densities derived in the following may 
or may not be related to density or column density th resholds for 
star format i on introduced in o bservational studies (e.g. , |Heiderman| 
|et al.||2010| |Lada et a.1. ||2010] ) . Studying such potential relations, 
nowever, deserves lurthcr consideration in the near future. 
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multi-freefall versions of the aforementioned theories, yet 
without considering magnetic fields. We build on all 
these approaches and extend the non-magnetic multi- 
freefall models to include magnetic fields. We then deter- 
mine the best combination of the aforementioned theo- 
retical ideas to come up with a more universal theoretical 
model for the SFR. Table Q] summarizes all six theoreti- 
cal models, which are discussed and derived in detail in 
the following. 

2.4.1. The KM Model 



In the KM model by |Krumholz fc McKee| j 2005[ ), the 
freefall-time factor tg ( po)/tff(p) in Equatio n ffil is simply 
set to unity. Moreover, Krumholz fc McKee| ( |2D05| ) define 
the critical density 5 cr j t i n the lower li mit of the SFRff 
integral by comparing the Jeans (1902) length 



Aj(p) 



(ncl 



\Gp 



1/2 



(9) 



evaluated at the me an density with the sonic scale A s 
(defined in Equation 13 below), 



^crit 



2 In 



(10) 



This choice is motivated by the expectation that the col- 
lapse sets in roughly at the sonic scale, where the tur- 
bulent fluctuations are of the order of the thermal sound 
speed, i.e., the local Mac h number has dropped to about 



unity at the sonic scal e (|Vazquez-Semadeni et al. 2003 
Federrath et al.|2010b| ). The global turbulent supersonic 
support is expected to become insignificant at the sonic 
scale, such that collapse can pr oceed below that scale 



Mac Low fc Klessen || 2004 |. The leading factor 2 
10| stems from the density dependence of 



(e.g. _ 
in Equation 

the Jeans length, Aj(p) oc p" 1 ^ 2 , and the numerical fac 
tor <j) x allows for slight vari ations in the actual s c ale on 
which the collapse sets in. Krumholz fc McKee|( 2005) 
find <j> T . = 1.12 for the simulations by |Vazquez-Semadeni| 



et al. (2003). In real molecular clouds, the sonic scale is 



Falgarone et al. 1992; Goodman et al. 


1998; Stahlcr & 


Palia|2004| 


Schnee et al. 2007;;|McKee & Ostriker|2007). 


To make 


Equation (|lo| more useful, we express all de- 



bers. This can be achieved by rewriting the Jeans length 
as 



Aj(po) 



= ( g£g 
\Gp 



1/2 



6GM r 



1/2 



(11) 



where we have assumed a spherical cloud with diameter 
L, mass M c , and isothermal sound speed c s . Since the 
velocity fluctuations in a turbulent medium depend on 
the length scale I as 



,(£)=a v (e/L)", 



(12) 



where ay ~ lkms -1 is the three-dimensional, non- 
thermal velocity dispersion on the scale L ss 1 pc, and 
p » 0.5 from observations in Galactic clouds ( Larson 



1981 ; Solomon et al. 


1987; Osscnkopf & Mac Low||2002 


Heyer & Brunt 200^ 


t Heyer et al. 


2009 Roman-Duval 


et al.|201ip, the Galactic Gentral JV 


Loiecular Zone <\ J ones 



et al.| |2012; Sh etty et al . [2012 ), and from numeric al sim- 



ulations dRritsiik e t al.||2007[ |Schmidt et al.|2009| |Feder- 
rath et al. 2010b), the sonic scale can be written as 



A s = L (Cs/fJy) 17 



(13) 



Substituting Equations ( 11 1 and ( 13 ) into Equation ( 10 1, 
we find 



Scrit.KM 



hi 



5 6GM C V c s 



5a v L ( ay 



2{l-p)/p 



= ln[(7r 2 /5)^, 



;M 2 



(14) 



where we have identified the virial parameter for a spheri- 
cal, uniform-density cloud with velocity dispersion ay on 
the diameter scale L, 



5a v L/(6GM c ) , 



(15) 



and the rms Mach number, A4 — ay / c s , and used p = 0.5 
in the second step. Thi s derivation is essentially ide ntical 
to the one presented in |Krumholz fc McKee| (1 2005), with 
the exception that we use the more general expression 
for the virial parameter here, 



2-Ekin/ I -Eg 



(16) 



the ratio of twice the kinetic energy to the gravita- 
tional energy. This general form reduces to a v ir,o given 



by Equation (|15| with E^ n = M c a v /2 and E gTS 



— 3GM c 2 /(5i?) for a spherical, homogeneous cloud with 
radius R = L/2. We emphasize that the definition of 
a V ir,o is based on global parameters, assuming a spheri- 
cal cloud with uniform density. This is far from realistic, 
given that clouds are in fact highly inhomogeneous and 
non-spherical. In the analytic de riva tions, however, this 
simplification given by Equation ( |15[ ) is necessary to en- 
able a mathematical analysis of the problem. In the sim- 
ulations discussed in Section [3] below, however, we will 
directly compute the virial parameter from the gravita- 
tional potential of the actual, three-dimensional, inhomo- 
geneous spatial gas distribution, providing a more gen- 
eral and accurate measure of the virial parameter given 
by the general form, Equation ( 16 ). This is discussed fur- 
ther below when we compare the theories to numerical 
simulations and in Section |7. 1.31 

The original model by 
gleets magnetic fields 



Krumholz fc McKee| ( |2005[ ) ne- 
iere, magnetic-field effects are 
partially added automatically by using Equation Q for 
a s , such that a s de creases with increasi ng magnetic en- 
ergy as derived in Molina et al. (2012). This however 
only changes a s , while a modification of s cr ;t is also nec- 
essary to fully account for magnetic-pressure effects on 
SFRff. 

Here we provide and apply a simple rule to include 
magnetic-field effects in the expression for the critical 
density s cr it . The key idea is to replace the thermal pres- 
sure by the sum of the thermal and magnetic pressures: 



2 



' Pth 
2 



Prr 



P c;^pc; + {l/2)pv A 



(17) 



where the second line implies isothermal gas. Using v\ = 
2c 2 /3 _1 with the definition of plasma f3 = Pth/-Pmag in 
Section 2.2 we can thus simply replace the sound speed 



Fcdcrrath & Klessen 



TABLE 1 

Six Analytic Models for the Star Formation Rate per Freefall Time. 



Analytic 
Model 


Freef all-time 
Factor 


Critical Density p C rit/po = 


= exp(s crit ) 


SFRg 






KM 


1 


(tt 2 /5)</> 2 xa vir M 2 (l + 




e/(2*t){l- 


t- erf [(a 2 - 


2s crit )/(8cr 3 2 )V 2 ]} 


PN 


*ff (Po)Aff(Pcrit) 


(0.067) 8- 2 xa vil M 2 f(/3) 




e/(2*t){l- 


(- erf [(a 2 — 


2 Scrit )/(8 CTs 2 ) 1 /2]}exp[(l/2)s crit ] 


HC 


tff(po)/tffO) 


(n 2 /5)y^ t xa viT M- 2 {l 


+ /3 1 ) + Pcrit.turb 


e/(2«t){l- 


^erf [(<x 2 - 


Scrit )/(2^) 1/2 ]}exp [ (3 /8)<x 2 ] 


multi-ff KM 


tff(po)/*ffO) 


(tt 2 /5)</> 2 xa vil M 2 (l + 




e/(2*t){l- 


t" erf [(°"? ~ 


Scrit )/(2 CT 2)V2]}exp [(3/8)<x 2 ] 


multi-ff PN 


tff(po)/*fi(p) 


(0.067) e~ 2 xa vil M 2 f(/3) 




e/(2«t){l- 


r erf [(<x 2 - 


Scrit )/(2 CT 2)V2]}exp [(3/8)<x 2 ] 


multi-ff HC 


t«{po)/tg(fi) 


(7r 2 /5)y: 2 t xa vir M- 2 (l 




e/(20 t ){l- 


Herf [(<x 2 - 


Scrit )/(2a 2 )i/2]} ex p [(3/8)<x s 2 ] 



Notes. The function f(fi), entering the critical density in the PN and multi-ff PN mod els is given by Equation ( |31| 
turbulent contribution p cr it,turb m the critical density of the HC model is given by Equation ||39l. 



The added 



by an effective sound speed, 

c s -> c s (f + 13 



-IN V2 



(18) 



Since M. = av/c s , we can also replace the sonic Mach 
number by an effective Mach number to take magnetic 
pressure into account: 



M ->M (l + /T 1 ) 



-1/2 



(19) 



Doing this for s cr it,KM in Equation (14 1 yields the mag- 
netic version of the critical density, 



s C rit,KM = ln (n 2 /5)(j) 2 x a vir M 2 (l + /3 1 )" 



(20) 



Even though we simply replaced the thermal sound speed 
by an effective, magnetic sound speed to derive this ex- 
pression, it has a deeper physical meaning. What we 
physically do in the derivation of s crrt is to replace the 
thermal Jeans length in the numerator of Equation ( 10 1 
with the magnetothermal Jeans length, 



* J, mag 



(1 + 0- 



1/2 



Gp 



(21) 



and the sonic scale in the denominator with the magne- 
tosonic scale, 



i/p 



(22) 



We note that the magnetic modifications given by Equa- 
tions (17) only account for magnetic pressure, i.e., 



isotropic pressure induced by the small-scale magnetic 
field. It does not account for mean magnetic-field effects, 
and as such will only be a valid extension to MHD as long 
as the turbulence remains trans- to super- Alfvenic be- 
cause sub-Alfvenic turbulence with a strong mean mag- 
netic field component is anisotropic, which is discussed 
at more detail below. 

Finally, solving the general SFRg-integral (Equation[7| 
with s cr it = Scrit, km from Equation ( 20 ) and unity for trie 
freefall-time factor (see Table [l] for a summary), the SFR 
per freefall time in the KM model is 



SFR, 



ff.KM 



exp(s)p(s) ds 



e 



1 + erf 



a 2 - 2s 



crit.KM 



(23) 



This deriva tion is identical to the one in lKrumholz fc Mc^l 



Kee| ( 2005[ ) , except for the extension to include magnetic 
fields in the theory based on the plasma j3 terms in <j s , 
Equation Q , and in the critical density, Equation ( 20 1 . 

2.4.2. The PN Model 

|Padoan fc Nordlund ( |2011 ) use ts(po)/tg (p C rit) as the 
freefall-time factor tg(po)/tff(p) in Equation Q, such 
that the freefall time of the critical density is used for 
all den sities above the critical dens ity to estimate SFRg . 
Unlike Krumholz & McKee (2005) who relate the criti- 
cal density s rr , t to the Jea ns length and the sonic scale, 
Padoan & Nordlund (2011 ) related the critical density to 
the magnetic shock jump conditions and to the magnetic 
critical mass for collapse. Starting with their assumed 
balance of thermal plus magnetic pressure by turbulent 
ram pressure on the cloud scale, 



/'Aim-) ( c s + \ v a) = Po ('y) 



(24) 



and using the definitions for Ai and (3 from Section 2.2 
Padoan & Nordlund (2011) arrive at an expression for 
tfie density jump 



Pmhd — Po - 



M 2 (3 
4~/3TT 



This leads to the post-shock thickness 

4 (3 + 1 



A 



MHD 



= OL 



M 2 p 



(25) 



(26) 



since Pmku/Po = ^-^/Amhd with the numerical param- 
eter 6 < 1, the fraction of the cloud size forming the 
largest shocks. Thus, 9L can be interpreted as the turbu- 
lent injection or forcing scale. In numerical simulations, 
most of the kinetic, turbulent energy is usually injected 
at a wavenumber k — 2 in units of 2tt/L, corresponding 



to half of the total cloud size 



Schmidt et al. 2009 Fcderrath et al. 



Kritsuk et al. 2007 



010b), as in the sim- 



ulations discussed below in Section|5] Thus, 9 « 1/2, but 
there might be some c o rrections to that particular scale 
fWang fc ; George|2002j) . |Padoan fc Nordlundj poTTj ) take 
v ~ 0.35. Here, we will simply interpret 9 as a numerical 
factor of order unity, accounting for any uncertainty in 
the post-shock thickness with respect to the total cloud 
scale L in Equation ( 26 ) . 



In order to derive a criti cal density for star formation, 
Padoan & Nordlund ( 2011 ) compare the mass of a sphere 
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with radiu s Amhd/2 to the critical mass for collapse. Mc- 
Kee| ( |1989[ ) define the critical mass for collapse of a mag- 
netized gas sphere as 



where 



Af crit « A/be + M$ , 
M BE = l.l82c 3 s G- 3/2 p- 1/2 



(27) 



(28) 



is the Bonnor-Ebert mass (Ebert 1955 Bonnor 1956 1 and 



nR 2 B 
G 1 / 2 



9tt 5 / 2 

'2G3/2 



p- 1/2 vi 



(29) 



is the magnetic critical mass for a sphere with radius R, 
threaded by a magnetic field B, where we have used the 
Alfven velocity «a = B /(inp) 1 / 2 in the second step. The 
numerical factor m$ in Equation ( 29 ) can vary depending 
on the geom etry and model taken, e.g., Pado an fc Nord- 
lund1(|20rfl) take rn$ = 0.17 with a reference to |Tomisaka| 
et al.f qi988D while |McKee| ( |1989[ ) use m<j> = 0.12, and 
Strittmatter| ( |l966[ ) derive = (127r 2 /5) _1/2 ~ 0.21 for 
a non-rotating cloud and m$ = (97r 4 /10) -1 / 2 ~ 0.11 for 
an oblate spheroidal cloud with eccent ricity approaching 
unity (see |Nakano fc Nakamura||1978 [). 
Finally, inserting Equations (25) and (29) into Equa- 

^Wcrit(Pcrit) = 



tion ( |27[ ) and setting the critical mass 
(47r/3)(AMHn/2) 3 p n rit, with the post-shock thi< 
given by Equation (|26[), yields the critical density, 



s C rit,PN = In [0.067e- 2 a vil M 2 f((3)] (30) 



with 



(i + 0.925/3- 3 / 2 
(1+r 1 ) 2 



2/3 



Note that s cr it,pN has the s ame dependence on a v 
M. as Scrit.KM in Equa t ion ( 20 ) . 



(31) 



and 



Padoan & Nordlund (2011 (use a rather special defini- 
tion of ji 1 which is the average post-shock /3. From a semi- 
analytical comparison of the mean magnetic field with 
the rms magnetic field, they derive a criterion for /? based 
on the av e rage A lfven Mach number, which Padoan fc| 
Nordlund (20111 simply use as a switch between MHD 
and purely HU turbulence. However, it is not straightfor- 
ward to derive a post-shock value of because it involves 



a density-threshold dependence (see discussion in |Padoan 
fc Nordlund||2011|). Moreo ver, the switch discussed by 



Padoan & Nordlund (2011) is a semi-analytical criterion 
derived trom their simulations. We therefore decide to 
ignore this spe cial definition of f3 for simplicity and apply 
Equation (30) with our definition of /3 (see Section 2.2), 
which includes all, and not just the post-shock gas. This 
is consistent with the definition of all other dynamical 
quantities of interest, e.g., a v ir, At, M.a, Po, etc. 

Using i ff (po)Aff (Pcrit.PN) and inserting s crit ,pN into the 
general Equation |7|) for SFRff yields 



SFR. ; .\ = y exp [ \s C]lA ,k\ 



exp(s) p(s) ds 



1 

— exp -S cr it,PN 

2(p t \2 



Scrit.PN 

2 



1 +erf 



PX 



(32) 



for the PN model. 



2.4.3. The HC Model 



Hennebelle & Chabrier (2011) were the first to argue 
that the treetall-time factor tff{p )/ts(p) must be used in 
Equation Q, such that different densities contribute to 
SFRff with their individual freefall time (see Equation[8]). 
The full HC model for SFRff is based on the mass spec- 
trum of gravitationally bound struc tures, as derived in 
Hennebelle & Chabrie'rl ((2C^1 [20091 : 



Af(M) = 



d(N/V) 
dM 



1 ds 
M~dM 



exp(s)p(s), (33) 



which is essen tially Equation (6) in Hennebelle fc| 



Chabrierl (1201 1|), except for the freefall time factor. The 
SFR in the~F 



HC model is then given by the integral over 
the mass spectrum, weighted by the mass and the freefall 
time factor: 



SFRff 



M dM ds tg(po) 
M dM t B (p) 



exp(s)p(s) 



tajpo) P 
tfi(p) Po 



p(s) ds . 



(34) 



N ote that the first equality is th e same as Equation (7) 

(2011_F It can be simplified 



Hennebelle & Chabrier 



to the second line in Equation (o4), by transforming the 



mass variable into the logarithmic density variable s and 
changing the limits of the integral accordingl y. W e em- 
phasize that the second equality in Equation (34) is ex- 



actly the same as the general model for SFRff given by 
Equation fff\ above. 

In the HC model, the critical density s C rit,HC is de- 
fined by requiring that the turbulent Jeans length Aj. tU rb 
at the cri tical density is a fraction y cut of the cloud 



scale L. Hennebelle & Chabrier (2011) do not pro 



vide an explicit physical interpretation of this choice, 
but a follow-up study is in preparation (P. Hennebelle 
& G. Chabrier 2012, private communication). The tur- 
bulent Jeans length is obtained by adding an effective 
turbulent press ure (see Chandrasekhar|[l951a|b[ Bonaz- 
zola et al.|l987jP1 to the sound speed in the purely thermal 
Jeans length, Equation 



M,turb 



_ ( 71-C 2 + (7r/3)gg(Aj ;t urb) 

Gp 

7tc 2 + 7rA Jiturb o- 2 / /(3L) 
G P 



1/2 



1/2 



(35) 



in which the turbulent velocity dispersion, Equation ( 12 ), 
must be evaluated on the scale t — Ajturb, such that 



5 Equation (7) in |Hennebelle &: Chabrier] |201l| | contains an er- 
ror in that the factor &1V1/M in tneir integral must instead read 
dAf (P. Hennebelle & G. Chabrier 2012, private communication), 
which simplifies the equation significantly because the mass and 
radius dependencies drop entirely and the integral can be com- 
pletely r ewri tten in terms of s and solved analytically (see our 
Equation [34]) . 

6 The concept of turbulent pressure is also used to derive ac- 
cretion rates and lumino sities during high-mass sta r formation in 

massive turbulent cores (IMcKee & Tan 2002, 20031. 
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the t urbu lent Jeans length is implicitly defined by Equa- 
tion ( 35 ) . Rewriting yields a quadratic equation with 
two solutions: 



A 



J.turb 



(P) = 



tt<4± y/3&irc%GL 2 p + 



QGLp 



(36) 



for which only the positive root is physical because the 
Jeans length must become larger when adding a stabi- 
lizing pressure — in this case a turbulent pressure. Natu- 
rally, this expression reduces to the thermal Jeans length 
for ay — » 0. Now, setting th e turbulent Jeans length 
equal to y cu tL as defined in |Hennebelle fc Chabrier| 



(2011), and identifying the virial parameter, Equa- 



tion 
yields' 



15), and the Mach number M. = oy/c s , finally 
The critical density threshold in the HC model: 



Scrit, HC = m [Pcrit,th + pcrit.turb] j 

where the (magneto)thermal contribution is 

p^th = (ir 2 /5)y~* t a vil M- 2 (l + /T 1 ) , 
and the turbulent contribution is 

/5crit,turb = (tT 2 /15) y~ u \ Q!vir • 



(37) 



(38) 



dr • (39) 

Note that the dependence of the thermal contribution 
to Scrit.HC on c^vir is the same as in the KM and PN 
models, while the dependence on the Mach number is 
Mr 2 , which is the opposite of the dependence in the 
KM and PN models, for both of which /o cr ; t oc M +2 (see 
Table [I] for a summary of all analytic models). 

The original HC model does not take magnetic fields 
into account, but we have extended the HC theory to 
MHD here by replacing the sonic Mach number in Equa- 
tion (38) with the magnetic version in t he same way as 
done tor the KM model via Equation (19). The mag- 
netic correction factor (1 + in Equation (38) simply 



becomes unity in the hydrodynamical limit [p — > oo) 

The SFR in the HC model is thus given by integrating 
Equation (34) or equivalently Equation ^ with s cr it,HCj 
which yields 



SFR ff . 



HC : 



S C rit,HC 

< /3 



ex P ( 2 s ) P( s ) ds 



20. 



exp —a 



1 + erf 



0~ s ~ s crit,HC 



(40) 



2.4.4. The Multi-freefall KM Model 



Following Hennebelle & Chabrier (2011), we define all 
three multi-lreefall versions of the KM, PN, and HC 
models by solving the generalized, multi-freefall integral, 
Equation ([7| . The analytic solution of that equation for 
an arbitrary threshold s cr i t is 



SFR ff = — exp 
2<?>t 



1 +erf 



> (41) 



which is identi cal to Equation (8) in Hennebelle & 
Chab rier (2011), and identical to the HC model, Equa- 
tion (3D), except that the critical density is defined ac- 
cording to either the KM, PN, or HC models. Thus, the 



multi-ff KM model is defined by using the threshold den- 
sity s cr it = Scrit .km from Equation (20) in the g ene ralized 
solution of the multi-freefall SFRff ,T?quation (41). 



2.4.5. The Multi-freefall PN Model 

The multi-ff PN model is defined by u sing the threshold 
density s crit = s cr ;t,pN from Equation (30) in the gener- 
alized solution of the multi-freefall SFRff, kquation (41 ). 



2.4.6. The Multi-freefall HC Model 

The multi-ff HC model is defined by taking the thresh- 
old density s cr it = s cr it,Hc from Equation (37), but only 
with the thermal contribution /5 C rit,th from Equation ( 38 ), 
while setting the turbulent contribution /5 C rit,turb = U, 
and using that threshold density in the ge nera lized so- 
lution of the multi-freefall SFRff, Equatio n j4lj). We do 
this to be consis tent with the definition in Hcimcb elle fcl 
Chabrier (2011). Note that the thermal density thresh- 
old is derived by requiring that the thermal Jeans length 
at that density, Aj is equal to y cn tL, while the full HC 
model includes the turbulent contribution, which is ob- 
tained by setting Aj_t U rb = ycutL (see the derivation of 
the HC model above). 

2.5. Dependencies of the Analytically Derived SFRff 

After the detailed derivation of the six different SFRff 
models, we can now start to compare them. Figure [l] 
shows all six SFR ff models: KM, PN, HC (left panelsjT 
and multi-ff KM, PN, HC (right panels) for a turbulent 
forcing parameter b = 0.4, corresponding to a statistical 
mixture of sole noidal and compressive modes in the tur- 
bulent forcing ( |Federrath et al.||2010bl Figure 8). When 
looking at the derivations above, it becomes clear that 
SFRff is a function of a v - lr , M, b, and f3. The depen- 
dencies enter through the definition of the critical den- 
sities in the different models, and through the variance 
of turbulent density fluctuations, Equation (|4). We plot 
the analytically derived SFRff as a function ra the virial 
parameter a v i r and the sonic Mach number M. in each 
panel. Note that all these models are plotted for (3 — » oo, 
i.e., without taking magnetic fields into account yet. As 
shown in Table [TJ each model has two fudge factors of 
order unity. The first one is l/c6 ( for all models (where 
the local efficiency was set to e = 1 for simplicity in 
all models, to facilitate the comparison), while the sec- 
ond one is <fi x , 6, and y cu t for the (multi-freefall) KM, 
PN, and HC models, respectively. We plot all curves 
for e/4>t — 1 to enable direct comparisons, and used the 
favored values of the fudge factors by the di fferent au- 
thors, (j) x = 1.12 j[R7u"mliolz fc Mc jKee]|2005 ]), 9 = 0.35 
Padoan fc Nordlund 2011), and y cu t = 0.1 (Hennebelle 
Chabrier||201l[ ). 

Dependence on a vir — Let us first concentrate on the de- 
pendence of SFRff on the viria l parameter. Since the 
virial parameter, Equations ( 15 1 and ( 16 ), is defined here 



as the ratio of twice the kinetic to the gravitational en- 
ergy, it essentially measures how strongly the system is 
bound, and whether it is contracting (a V ir ^ 1) or ex- 
panding (a V j r > 1). Thus, we generally expect that the 
SFR should decrease with increasing a v i r because the 
cloud then becomes less bound and less likely to form 
stars. Indeed, the analytic SFR generally decreases with 
increasing a v ir in all models with the exception of the 



The Star Formation Rate 9 



KM ($3=1.12) multi-ff KM ($ x =1 .1 2) 




Fig. 1. — Comparison of the six analytic models for the star formation rate per fr eefal l time, SFRff : KM, PN, HC (left panels) and 
multi-freefall KM, PN, HC (right panels). See Tableland the derivations in Section |2.3| for details of the different analytic models and 
functions plotted (e/(f>t = 1 in each panel). The dependence of SFRff on the virial parameter a v f r , and the sonic Mach number M are 
shown in each panel for a turbulent forcing parameter b = 0.4, corresponding to a statistical mixture of solenoidal and compressive modes 
in the turbulent forcing. All models are plotted without taking magnetic fields into account, i.e., plasma /3 — > oo. 
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original PN model, for which SFRff first increases for 
«vir ^ 1 and then decreases for a v i r ^ 1- The increase 
comes from the freefall-time factor iff (po)Aff(Pcrit) in 
the PN model, which leads to the factor exp(s cr it,PN/2) 
in Equati on (32), and with the critical density from 
Equation (30) to SFRff oc a v ; r for small ot v - K . As ex- 
pected though, this direct proportionality disappears in 
the multi-freefall PN model, as in the other two multi- 
freefall models (multi-ff KM and multi-ff HC). 

Dependence on M — The expected dependence of SFRff 
on the sonic Mach number is that SFRff should increase 
with increasing M. because higher Mach number means 
stronger and denser local compression, leading to higher 
SFRs. Indeed, the Mach number dependence is gener- 
ally similar in all models, i.e., SFRff increases with A4, 
with the exception of the original KM model, which has 
the weakest dependence on M.. For large a v i r , SFRff km 
increases, but only slowly, while for small a v i r , it stays 
constant or even decreases with increasing M. . Both the 
HC and multi-freefall HC models have the strongest pos- 
itive correlation with the Mach number, such that for 
M > 10, SFRff, hc hardly depend s on a vir anymore (see 
also, Hennebelle fc Chabrier|20TT| Figure 10 The strong 
increase of SFRff with M. in the two HC cases comes from 
the Mach number dependence of the thermal contribu- 
tion to s cr it,HC! which is p C rit,th Mr 2 , leading to a 
decreasing threshold density in the HC models, and thus 
to a higher SFR. This is the opposite compared to the 
KM and PN models, for both of which the critical den- 
sity increases with the square of the Mach number (see 
Equations 20 and 30 respectively, or Table [T]). 

We also note theTocal minima of SFRff around M ~ 2 
in all models, except the HC and multi-freefall HC mod- 
els. Those minima are spurious because they occur close 
to M — 1, for which the basic approach of shock-induced 
star formation must eventually break down as the system 
becomes transonic. Shocks require M > 1, by definition, 
but for rms Mach 1-2, a significant fraction of the sys- 
tem is transonic to subsonic. We thus conclude that all 
six models break down for the low Mach number regime, 
M, <2. The rms sonic Mach num bers in real molecular 
clouds usually exceed unity by far flLarson||1981[ |Falgar 
one et al.|1992 Roman-Duval et al.|2011||Schneider et a' 



2 F 



2012p , such that our analytic models are generally appli- 
cable to typical molecular clouds with M > 2. 

Dependence on b — While the dependence on Mach num- 
ber enters SFRff both through s cr ;t and a 2 , the forcing 
dependence only enters through the forcing parameter b 
in a 2 , Equation Figure [2jshows SFRff as a function 
of the forcing parameter b for all models and three differ- 
ent Mach numbers (M — 5, 10, and 20). All curves are 
plotted for a v - lr — 1, j3 — > go, e/(j) t — 1, and the standard 
fudge factors 4> x = 1.12, 9 = 0.35, and y cut = 0.1, respec- 
tively. We see that SFRff increases monotonically with 
b, from 6=1/3 (solenoidal forcing), over b = 0.4 (mixed 
forcing), to b = 1 (compressive forcing) in all models. 
This is expected because the density variance becomes 
larger for more compressive forcing, pushing a significant 

7 Note t hat the three different sonic M ach numbers shown in 
Figure 1 of [Hennebell e &: Chabrier| | |2011| are actually Ai = 4.5, 
9, and 18, and not 4, y, and Id as indicated in their figure caption 
(P. Hennebelle & G. Chabrier 2012, private communication). 
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Fig. 2. — SFRff as a function of the forcing parameter b in Equa- 
tion Q for sonic Mach numbers Jvi = 5 (top), M = 10 (mid- 
dle), and M = 20 (bottom). All curves are plotted for ct v i r = 1, 
t/d>t = 1, and the favored fudge fa ctors by Krumholz & McKee 
|Padoan fc Nordlund| ( |2011[ ) , and |Hennepl5eT^^naEfier 
<px = 1.12, y = 0.35, and j/ C ut = 0.1, respectively. Unly 
hydrodynamic cases are shown (/3 — > oo). The star forma- 
tion rate increases monotonically from 6=1/3 (solenoidal turbu- 
lent forcing), over b = 0.4 (mixed forcing), to b = 1 (compressive 
forcing) . 



fraction of the gas to higher densities 


Federrath et al. 


2008b 


2010b 




Konstandin et al. 


2012a 


. Similar to the 



behavior with increasing Mach number, increasing the 
amount of direct compression induced by the turbulent 
forcing leads to higher local densities, and thus to higher 
SFRs with a typical increase of about an order of mag- 
nitude for compressive forcing compared to solenoidal 
forcing. 

Dependence on fi — We expect that by adding magnetic 
energy to the system, the SFR should decrease because 
magnetic energy adds a stabilizing pressure to the sys- 
tem, counteracting gravitational collapse. Figure[3]shows 
the dependence of SFRff on plasma j3 in the six analytic 
models. We emphasize that only the original PN model 
had a magnetic-field dependence, co ming from the de- 
pendence of s cr it,pN on (3 in Equation (30), and from the 
dependence of a s on (3 in Equation (H). However, we have 
extended all other analytic models (KM, HC, and multi- 
ff KM, PN, HC) to MHD, simply by applying the MHD 
version of a s , Equation (|4| in all models, and replacing 
the sonic Mach number m the expressions for the criti- 
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Fig. 3. — Same as Figure [2] but SFR ff is shown as a function 
of plasma /3, the ratio of thermal to magnetic pressure (bottom 
abscissa) or as a function of the Alfven Mach number, Ma = 
M^J 13/2 (top abscissa) for mixed forcing (b = 0.4). Since the 
sonic Mach number is M = 5, 10, and 20 (top to bottom panels), 
the A^A-axis varies between the three panels. The solid, vertical 
line separates Ai\ < 2 fr om M a > 2. Analyt ic predictions below 
A4a $ 2 are inaccurate ( |Molina et al.|[2012[ l and only shown in 
gray. 



cal density by the magnetic version M. 
introduced in Equation (19 1. 



As found in a detailed comparison of the analytically- 
derived o~ s with numer i cal si mulations of MHD turbu- 
lence in 

breaks down for 



Molina et al. (|2012|), the jtandard deviation- 
Mach number relation, 



Equation (14 
2 because stron gly sub-AIfvenic flows become 
highly anisotropic (e.g., |Mac Low||1999[ |Cho & Vish- 



M 



< 



Brunt et a 



macl[2000l [Cho fc Lazarian||ii()03l |Beresnya.k et, al.P 0051 



2010| |Esguivel fc Lazarian||201ip . Since the 



magnetic-held dependence ot SKKff was introduced as 
an isotropic magnetic-pressure extension, the behavior of 



the analytic models for A^a < 2 is likely invalid. Thus, 
we only consider the trans- to super- Alfvenic regime with 
A^a ^ 2. In this regime, SFRff decreases with increasing 
magnetic energy, i.e., decreasing /3 or A^a in all models, 
as expected when adding a stabilizing magnetic pressure. 

3. TESTING THE ANALYTIC THEORIES FOR THE SFR 
WITH NUMERICAL SIMULATIONS 

In order to test the analytic predictions of the star for- 
mation rate (SFR) models in Section [2] we perform a 
series of numerical simulations of driven, supersonic tur- 
bulence, including magnetic fields, gravity, and a model 
for collapse and accretion of star-forming regions to mea- 
sure the SFR. Ideally, we would like to sample as much 
of the parameter space as possible with the simulations. 
Since t he an alytic SFR depends on a v i r , AJ, b, and /3 (see 
Section 2.5), we have to restrict ourselves to testing only 
a subset ot those because the simulations are computa- 
tionally too expensive to scan through the entire parame- 
ter range. We thus concentrate here on the Mach number 
and forcing dependence, as well as the dependence on the 
magnetic field, but only consider models with an initial 
virial parameter of around unity. However, as the tur- 
bulence produces strong spatial density variations, the 
virial parameter can change by an order of magnitude 
from its initial value given by Equation (15 1 when the 



turbulence is fully established because the mass is re- 
arranged into complex filamentary and sheet-like struc- 
tures. To take this into account, we always compute 
instantaneous values of a v \ ri ba sed on the spatial distri- 
bution of the gas (Equation 16), as for all other param- 
eters, and then average them over space and time. The 
time interval for averaging is chosen such that it cov- 
ers the whole star formation sequence in the simulations, 
from the time when the turbulence is fully established, 
as explained in more detail in Section |3.4| below. First 
however, we explain our numerical scheme in Section |3.1| 
the forcing of the turbulence in Section [3?2} and the sink 
particles introduced to model core and star formation in 
Section 13.31 

3.1. Numerical Methods 



We use the adaptive mesh r efinement (AMR, Berger fc 



| Colella||l989) code FLASH^] ( |Fryxell et aL||2000l |Dubey 
|et al.||2(J(J8[ ) in version 2.5 to integrate the ideal, three 
dimensional, MHD equations, including self-gravity, 



dp 
di 



V-(pv)=0, 



0_ 
Of 



\" I v (B 4 J )B - VP, + p (g + F stir ) , 



dE 
~dt 
dB 

~dt 



+ V- 



(E + P*) v 



= V x (v x B) 



(B • v) B 

47T 

VB = 0, 



pv • (g + F stir ) 



(42) 



where the gravitational acceleration of the gas g, is the 
sum of the self-gravity of the gas and the contribution of 
sink particles (a subgrid model for collapse and accretion 
of star-forming regions in the simulations, explained in 



8 http : //flash. uchicago . edu/site/f lashcode/ 
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Section 3.3 below): 

V z $ gas = AirGp 



g = -VfPgas + gsinks ■ 
1* 



(43) 



In the ideal MHD Equations (42 1, p, v, P* = P t 



th 



1/(8tt) |B|*, B, and E = pe int + (p/2) |v| z + 1/(8tt) |B| z 
denote gas density, velocity, pressure (thermal plus mag- 
netic), magnetic field, and total energy density (inter- 
nal plus kinetic, plus magnetic), respectively. The MHD 
equations are closed with a polytropic equation of state, 
Ph = c 2 p r with r — 1, such that the gas remains isother- 
mal with a constant sound speed c s = 0.2 km s -1 , corre- 
sponding to a temperature of T sa 11 K for gas with a 
mean molecular weight of 2.3. This is a reasonable ap- 



ity, over a wide range of densities (Wolfire et al. 1995 


Omukai et al. 


2005; Pavlovski et al. 


2006; 


Glover & 


Mac 
Hen- 


Low||2007a|b 


Glover et al.||2010 




ill et al. I20TTT 


nemann et a 


2012). Moreover, 

Vrm — ' i i • 


Gl 


over & Glark| (|2012|) 



ity. Reducing the metallicity of the gas by two orders of 
magnitude reduces the time-averaged SFR by less than a 
factor of two. Thus, our conclusions remain intact, even 
though we neglect the detailed chemistry, cooling and 
heating processes in molecular clouds in th is study. 

We solve the MHD Equations (42) on three- 
dimensional, periodic grids with maximum resolutions 
of 7V r 3 s = 128 3 -1024 3 grid points. These are all uniform- 
grid simulations, except for the N IOS — 1024 simulation, 
where we use a root grid with 512 3 cells and one level of 
AMR with a refinement criterion to ensure that the local 
Jeans lengths is covered with at least 32 grid cells, in or- 
der to resolve turbulent vorticit y and magnetic-field am- 
plification on the Jeans scale ( Sur et aL]|2010| |Federrath| 
et al .|2011c Turk et al. 1 2012). We use a positive-definite^ 
J Rie mann solver flBouchut et al.||2007| |2010[ |Waa-| 
2009|, which has been tested to r effici ency, robust- 



MH 



gan 



ness, and accuracy in |Waagan et al. (2011). This study 
shows that the MHD scheme keeps V-B errors at a negli 



gible level, and allows us to model extremely high-Mach 
turbulence without producing unphysical states. This is 
particularly important for this study because we model 
supersonic turbulence on the largest scales of molecular 
clouds with rms Mach numbers as high as M. ~ 50 and 
compressive forcing, which produces density contrasts by 
several orders of magnitude, sometimes between two ad- 
jacent grid cells because of multiple interactions of shocks 
and strong rarefaction waves, even before gravitational 
collapse sets in. Grid-based HD solvers often produce 
negative densities in such situations because of numer- 
ical post-shock oscillations. Such unphysical states are 
avoided by construct ion in the HLL3R Riemann scheme 
(Waagan et al. 2011) used here. The self-gravity of the 
gas, i.e., tfie gas-gas gravitational interaction (Equa- 



tion 



43 ) is computed using a multi -grid Poisson solver 



(the ELASH2.5 version discussed in Rickcr 2008), while 



the sink particle interactions are comput ed by direct N- 
body summation, as explained in S ect ion [3 . 3 1 below . We 
note that the gravitational potential $g as is computed 
with respect to the periodic boundary conditions speci- 
fied in the simulations. 



The ideal MHD Equations ( 42 ) do not contain any ex- 



However, any numerical scheme has an effective numeri- 
cal viscosity v and magnetic resistivity r) due to the neces- 
sary discretization of the MHD equations. Even though 
the numerical viscosity depends on the specifications of 
the algorithm, it can be used to mimic the effe cts of ex- 
plicit viscosity and resistivity ( Benzi et al. 2008[ ). It is im- 
portant to realize, though, that the kinematic and mag- 
netic Reynolds numbers that can be achieved wit h ideal 
MHD depend on t he grid resolution. As shown in |Feder-~ 



rath et al. (2011b), compressible, ideal MHD turbulence 
resolved with 128 J grid cells reaches kinematic Reynolds 
numbers Re = Lay/v « 1500 an d magnet i c Rey nolds 



Burgers (11948b seal- 



numbers Rm = Lay jr\ « 3000. For 
ing of the turbulence a v (£) cx l 1 / 2 (Equation |l2| with 

3/2 — 

p = 1/2) the Reynolds num bers scale cx Afro's as op- 
posed to Kolmogorov (1941) scaling of the turbulence, 
<J V (£) cx Z 11 ' 6 (Equation 12 with p = 1/3), leading to 

a Reynolds-number scaling oc N^Is- Thus, even in our 
highest resolution simulation with N rcs = 1024, we only 
achieve Reynolds numbers, Re ps (2.4-3.4) x 10 4 and 
Rm ps (4.8-6.8) x 10 4 , depending on the scaling of the 
turbulence. In summary, although the flows we model 
exhibit fully developed turbulence (Frisch 1995), their 
Reynolds number are still considerably smaller than the 
ones inferre d for real molecular clouds (see, e.g., Schober 
et al. 2012). We will thus study the resolution depen- 



dence of our results for the SFR below. 

3.2. Turbulent Forcing 

Previous numerical studies of non-driven turbulence 
have shown that supersonic turbulence decays in about 
a crossing time, irre spective of whether magnetic fields 



are included or not (Scalo fe Pumphrey 1982; Mac Low 
et~S1[l998| |Stone et al.| [19981 l M ac Low| " |T9M| . The 



observed presence of turbulence has thus lead to the 
conclusion that interstellar turbulence should be driven 
by some physical stirring mechanisms. Those mecha- 
nisms include supernova explosions and expanding, ion- 
izing shells from previous cycles of star formation (|Mc- 
| KeeJl989j|Krumholz et al.|2006[|Balsara et al.|2004l|Bre^ 
litschwerdt et al.|2009| |Peters et al.|2011||Goldbaum et aT 



|2011||Lee et al.||2012 1, gravitational collapse an d accre- 



tion of material! Vazquez-Semadeni et al.||l998[ |Klessen 
fe Hennebe llc 2010; Elmcgrccn & Burkcrt 2010; Vazquez 
Federrath et al.||2011c| |Robcrtson 
and galactic spiral-arm com pression 



.Semadeni et al.||20 10 
& Goldreich||2012| 



of HI clouds ( |Dob'bs fe Bonn cll 2008; |pobbs et al.|2008L 
and magnet orot at ional instability (MRl) (Piontek fe (Js- 

|. On smaller scales, jets 



triker|2007||Tamburro et al.|2009| 

and outflows from young s tellar objects have been sug- 



gested to drive turbulence (|Norman fe Si lk 1980; Baner 
ee et al.l|2007HNak amura & Li 2008 ; Cunningham et a 



plicit kinematic viscosity and magnetic resistivity terms 



2009| lUarroll et al |2010| [Wang et al.||2010[ ). Turbulence 
in high-redshift galaxies is also likely driven, by feedbac k 
from previous cycles of star formation ( |Green et al.|2010"| . 
A summary and comparison of driving mechanisms for 
interst ellar turbulence is pro v ided in|Mac Low fc Klessen| 
pOfjl] and |Elmegreen) |2009| . |Mac Low fc Klessen| ( |2004T 
conclude that expanding shells are likely the dominant 
driver of interstellar turbulen ce in the star-for ming parts 
of the Galaxy. More recently, Lee et al. (20121 also noted 
that the kinetic energy injected per unit time by star- 
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forming complexes via expansion of bubbles is about 2/3 
of the luminosity required to maintain the observed ve- 
locity dispersions, supporting the view that expanding 
bubbles driven by massive star clusters from previous 
star formation are a major driver of turbulence in the 
Milky Wa y (see e.g., the Cygnus X giant molecular cloud 
studied in | Schneider et aLj2011|). 



In index notation, this tensor reads 



kj k a 



It is important to realize that all these potential drivers 
(maybe with the exception of the MRI) are expected to 
primarily drive compressible modes in the velocity field, 
but do not directly excite solenoidal modes. However, 
even though the turbulence in molecular clouds might be 
driven compressively, solenoidal modes are indirectly ex- 
cited b y nonlinear interactions of multiple colliding shock 
fronts (|Vishniac| |1994[ [Sun fc Takayama||2003| |Kritsuk 
et al.||20U7{ ~|Fedcrrath et al. 2010b ) , Tryljaroclimt y, ro 



V% = C + (1 - C) V% = ( Sv + (1 - 2C) ^f- , (44) 

where Sij is the Kronecker symbol, and V^j = 5ij — 

kikj/k 2 and v\j = kikj/k 2 are the solenoidal and com- 
pressive projection operators, respectively. The ratio of 
compressive po wer to total power in F st i r can be derived 
from Equation ( 44 1 by evaluating the norm of the com- 



pressive component of the projection tensor and dividing 
it by the total injected power, resulting in 



(W) 2 



F u 



1 - 2C + 3C 2 ' 



(45) 



tation and s hear (|Del Sordo fe Brandenburg|[20TT|), and for three-dimension al space (iSchmidt et al.||2009| IPed 



2011b I, such that supersonic turbulence driven by even 



^^^h^^^^^^h ^^^^^tf^^^^^^^^^^^^^M^^^^M 1 V.' » V 111 V. V V 11111 VLlkJIV i»ua X"^ "IVV \ h_/ V111111V* V V K> >.A, A ■ U V/ \/ I.' 1 V V 1 

by visc osity ( |Mee fe Brandenburg||2006| |Federrath et al.| lerrath et al.||2010b|). The projection operator servesTo 



purely compressive forcing contains about half of its ki 
netic power in solenoidal modes and t he other half in 



compr essible modes in the inertial range (Federrath et al. 
[2010b | Figure 14) 



Modeling physical turbulent stirring mechanisms in nu- 
merical simulations requires assumptions about the spa- 
tial and temporal correlation of the turbulent forcing 
events. It is also still a matter of debate which of the 
physical mechanisms dominates the injection of turbu- 
lent energy on different cloud scales. Given these un- 
certainties, instead of trying to mimic one or more of 
the potential physical drivers of turbulence, we here use 
simulations of the so-called 'driven turbulence in a box'. 
From these simplified and idealized simulations, we can 
draw statistical conclusions about the role of turbulence 
for star formation, given average properties of a cloud 
(a v i r , M, b, and /3). In particular, our turbulent forcing 
approach allows us to evaluate the role of the mixture of 
velocity modes excited by a physical driver. 

In practice, the stochastic forcing term F st i r is applied 
as a source term in Equations ( 42 ) to drive turbulence 
in the simulations. F st i r only contains large-scale modes, 
1 < k < 3, where most of the power is injected at the 
k = 2 mode in Fourier space, which corresponds to half 
of the box size L in physical space. We thus model tur- 
bulent forcing on large s cales, as favored by molecular 



cloud observations 



Heyer et al. 2006 ; Brunt et al 
Roman-Duval et 



Ossenkopf fc Mac Low 2002 
TTacnsler et al. 2011 



2009 



al. 



201ip . Smaller scales, k > 3 are 
not affected directly by the forcing, such that turbulence 
can develop self-consistently on these scales. We use the 
Ornstein-Uhlenbeck (OU) process to model F st i r , which 
is a well-defined stoc hastic process with a finite auto 



et al 
force 



2006), 
neid in 



correlation timescale (Eswaran & Pope 1988 Schmidt 



leading to a smoothly varying stochastic 



m space and time. Details about the OU pro- 
ce ss and the forcing app li ed in this study can be found 



in ISchmidt et al.| (2009), Federrath et al. (2010b), and 



Konstandin et al.| (2012a). However, the essential point 
of our forcing approach is that we can adjust the mix- 
ture of solenoidal and compressive modes of F st j r . This 
is achieved by decomposing a given vector field with ran- 
dom mixtures into its solenoidal and compressive parts, 
by applying the projection tensor V_ ^(k) in Fourier space. 



construct a purely solenoidal force field by setting ( = 1, 
while for £ = 0, a purely compressive force field is ob- 
tained. Any combination of solenoidal and compressive 
modes can be constructed by choosing £ £ [0, 1]. Here 
we compare simulations with C — 1 ( s °l)i C = 1/2 (mix), 
and C = (comp). A detailed study of the forcing de- 
pendence of the ^-parameter entering the expression for 
the variance o f the density PDF, Equa tions Q and ([5]), 
is provided in Federrath et al. (2010b Figure 8), where 
they measure o as a function of the forcing parameter £. 

3.3. Sink Particles and Resolution Criteria 

In order to model collapse and accretion of star- 
forming gas in the simulations, we use a subgrid model 
called 'sin k particles', which is a method originally in- 
vented by Bate et al. ( 1995[ ) for Smoothed Particle Hy- 
drodynami cs, ancTT irst adopt ed for Eul erian, AMR sim- 
ulation s by Krumholz et al. (2004 1. In Krumholz et al. 
(20041, a Lagrangian sink particle is introduced, if the 



gas reaches a given density. However, sink particles are 
supposed to represent bound objects that are going into 
collapse, and thus, a density threshold as th e only cri- 
terion for si nk particle creation is insufficient dFederrath 
et al.|2010a[). Based on the ideas of |Bate et al.| (1995 ) and 
Krumholz et al.| ( |2004[) , we use an advanced AlVlK- based 
approach tor sink particles, in which only bound and col- 
lapsing gas is accreted, thus avoiding the creat ion of spu- 
rious sink pa rticles (for a detailed analysis, see 
et al.pOlOa l. The key feature of this approac 



Federrath 
r is to de 

fine a control volume around cells that exceed the density 
threshold set by th e resolution criterion to avoid artifi- 
cial fragmentation. Truelove et al. (19971 found that the 
Jeans length must be resolved with at least 4 grid cells 
to avoid artificial fragmentation, leading to a resolution- 
dependent density threshold criterion for the creation of 
sink particles: 

TTC 2 

Psink = 777^1 ) (46) 



4G; 



sink 



where the sink particle accretion radius r S j n k is set to 
2.5 grid-cell lengths at the maximum level of refinement, 
corresponding to half a Jeans length at /c s ink, such that 
the Jeans length is still resolved with 5 grid cells prior 
to potential sink particle creation to avoid artificial frag- 
mentation. Grid cells exceeding the density threshold 
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given by Equation (46), however, do not form sink parti- 
cles right away. First, a spherical control volume with ra- 
dius r s i n k is defined around the cell exceeding p s i n k within 
which additional checks for gravitational instability and 
collapse are performed. We check whether the gas 

• is on the highest level of refinement, 

• is converging from all directions in the rest frame 
of the central cell (negative radial velocity), 

• is at a local gravitational potential minimum, 



• is bound (|-E grav | > 

• is Jeans-unstable, and 



E, 



ih 



• is not within r s i n k of an existing sink particle. 
If all these checks are passed, a sink p article is created 



in the center of the control volume (see |Federrath et al. 
2010a). This procedure avoids spurious sink particle for 



mation, and allows us to trace only truly collapsing and 
star-forming gas. Given the checks above, it is clear that 
in some cases, a sink particle is not necessarily formed 
even though the density threshold is exceeded. This does 
not mean, however, that such gas would be subject to ar- 
tificial gravitational fragmentation. Since the checks did 
not allow sink particle creation, the gas in the control 
volume was not collapsing and/or not bound, so there is 
no need to worry about artificial fragmentation at this 
stage, even though the density threshold was exceeded. 
This can happen quite frequently in supersonic turbu- 
lence because shocks can push the gas density above the 
threshold, even though this gas is not necessarily gravi- 
tationally bound after the shock passage. 

Once a sink particle is created, it can gain mass by 
accreting gas from the AMR grid, but only if this gas 
exceeds the threshold density, is inside the sink particle 
accretion radius, is bound to the particle, and is collaps- 
ing toward it. If all these criteria are fulfilled, the ex- 
cess mass above the density threshold defined by Equa- 
tion ( 46 ) is removed from the MHD system and added to 
the sink particle, such that mass, momentum a nd angular 
momentum are conserved by construction (see |Federrath| 
et al.||2010a||2011a| for details). 

All contributions to the gravitational interactions be- 
tween the gas on the grid and the sink particles are com- 
puted by direct iV-body summation over all grid cells and 
sink particles (gas-sink, sink-gas, and sink-sink), using 
gravitational spline softening inside the sink particle ra- 
dius to avoid singularities during close encounters. The 
softening only affects scales that are anyway below the 
grid-resolution cutoff set by the sink particle accretion 
radius. A second-order accurate Leapfrog integrator is 
used to advance the sink particles on a time step that 
allows us to resolve close and highly eccentric orbits of 
sink particles without introducing significant errors on 
super-resolution grid scales. 

3.4. Initial Conditions, Procedures, and List of Models 
Starting from a uniform density distributi on a nd zero 



velocities, the forcing term F st i r in Equations (42) excites 
turbulent motions. First, we evolve the MHD equations 
for two turbulent crossing times, 2T = L/{M.c s ) without 



pressible turbulence (e.g., 


Klessen et al.[ 2000 


Klessen 


2001; Heitsch et al.||2001l J 


Kederrath et al.||2009 


2010b; 


Price & Federrath||2010| \M 


icic et al. 


2012}. We do not 



that our measurements of the SFR are contaminated by 
this rather artificial initial transient phase, duri ng which 
the system is building up a turbulent cascade ( Schmidt 
et al.| 2009[ ). Af ter that, we solve the full system of MHD 



Equations ( 42 ) and ( 43 ) including self-gravity and for- 
mation of sink particles. For practical purposes, we reset 
the time t = 2 T to t = iff (po), which is the time when 
turbulence is fully established and star formation is al- 
lowed to proceed. We note that this procedure is slightly 
different from setting up a simulation with power-law ve- 
locity scaling drawn from Gaussian random seeds as an 



mation studies (e.g.,|Bate et al. |2003 


Clark et al.||2005 


Krumholz et al.||2007| |Frice & Bate 


2008, 2009 


Smith 


et al. 20081 Federrath et al. 2010a 


Walcn et al. 


2010 


Girichidis et al. 20111). In those cases, the initic 


d ran- 



profile (often constant density or radial power-law dis- 
tributions) , such that density and velocity fields have no 
causal connection. Here, the initial density and velocity 
fields at t = are consistently coupled via the equations 
of (magneto)hydrodynamics. We also keep driving the 
turbulence instead of imposing only an initial Gaussian 
perturbation as in the studies mentioned above. 

All our numerical simulations and their basic param- 
eters are listed in Table [2] Each model has a unique 
name, starting with 'GT' (for 'GravTurb'), followed by 
the maximum grid resolution ('128', '256', '512', and 
'1024'), the forcing type ('s':solenoidal, 'm':mixed, and 
'^compressive), and the Mach number ('M3', 'M5', 
'M10', 'M20', and 'M50'). Models with an initially uni- 
form magnetic field in the z-direction through the simu- 
lation box are additionally denoted with 'BE, 'B3', and 
'B10', corresponding to B — 1, 3, and 10/iG, respec- 
tively. Different random sequences with the same statis- 
tical properties for the turbulent forcing are indicated by 
'(si)', '(s2)', and '(s3)' at the end of the model name, in- 
dicating that random '(seedl)', '(seed2)', or '(seed3)' was 
used. Columns 2-10 in Table[2]list the maximum numer- 
ical resolution, type of forcing, mean density po, box size 
L, the total mass M c , large-scale velocity dispersion cry, 
initial magnetic-field strength Bq, initial plasma /3p, and 
virial parameter a V i r:0 computed with Equation (15). 

Columns 11-15 are derived quantities, measured as 
space and time averages after turbulence is fully estab- 
lished, t > 0, until 20% of the original cloud mass is 
accreted onto sink particles, i.e., the star formation effi- 
ciency has reached SFE = 20%. We list the average virial 
parameter a v - u -, the sonic Mach number M, forcing pa- 
rameter b, plasma (5, and Alfven Mach number Ma- The 
instantaneous virial parameter, Equation (16), in column 
11 of Table [2] is computed as a vu — 2E^ n /\E srav \ = 
^2 MivH\ Mi®gas,i\ from the gravitational pot ential 
$ ffa , a returned by the Poisson solver (see Section 3.1 1, as a 



sum over all grid cells i with mass M$ and velocity Vi. We 
note that this is dif ferent from the value a v ir,o obtained 
from Equation ( 15 ) and listed in column 10, which as- 
sumes a homogenous, spherical density distribution. In 
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TABLE 2 

Basic Parameters of the Numerical Models of Forced, Supersonic, Self-gravitating, (M)HD Turbulence. 



Model 




iVres 


Forcing 


P0 




L 


M c 


cry 


Bo 




QVir,o 


Qvir 


M 


b 




Ma 












[g cm - 


3 ] 


[PC] 


Wq] 


[kms- 1 ] 


[fiG] 
















(1) 






(2) 


(3) 


(4) 




(5) 


(6) 


(7) 


(8) 


(9) 


(10) 


(11) 


(12) 


(13) 


(14) 


(15) 


01) 


GT256sM3 




256 


sol 


5.8x10 


-19 


3.3X10" 1 


3.1 xlO 2 


0.59 





CO 


0.07 


1.4 


2.9 


1/3 


oo 


CO 


02) 


GT512sM3 




512 


sol 


5.8x10 


-19 


3.3X10" 1 


3.1 xlO 2 


0.59 





oo 


0.07 


1.4 


3.0 


1/3 


CO 


oo 


03) 


GT256mM3 




256 


mix 


5.8x10 


-19 


3.3X10- 1 


3.1 xlO 2 


0.61 





CO 


0.08 


1.1 


3.1 


0,1 


OO 


oo 


04) 


GT256cM3 




256 


comp 


5.8x10 


-19 


3.3X10" 1 


3.1 xlO 2 


0.58 





oo 


0.07 


0.46 


2.9 


1 


oo 


oo 


05) 


GT512cM3 




512 


comp 


5.8x10 


-111 


3.3X10" 1 


3.1 xlO 2 


0.58 





oo 


0.07 


0.48 


2.9 


1 


oo 


oo 


06) 


GT256sM5 




256 


sol 


3.3x10 


-21 


2.0x10° 


3.9xl0 2 


1.0 





oo 


1.0 


8.0 


5.0 


1/3 


oo 


oo 


07) 


GT256mM5 




256 


mix 


3.3x10 


-21 


2.0x10° 


3.9xl0 2 


0.99 





CO 


0.98 


5.4 


5.0 


0.4 


oo 


oo 


08) 


GT256cM5 




256 


comp 


3.3x10 


-21 


2.0x10° 


3.9 XlO 2 


0.91 





OO 


0.82 


1.5 


4.5 


1 


oo 


CO 


09) 


GT128sM10 




128 


sol 


8.2x10 


-22 


8.0x10° 


6.2 xlO 3 


2.1 





oo 


1.1 


11. 


10. 


1/3 


oo 


CO 


10) 


GT256sM10 




256 


sol 


8.2x10 


-22 


8.0x10° 


6.2 xlO 3 


2.1 





oo 


1.1 


12. 


10. 


1/3 


oo 


oo 


11) 


GT512sM10 




512 


sol 


8.2x10 


-22 


8.0x10° 


6.2 xlO 3 


2.1 





oo 


1.1 


12. 


10. 


1/3 


oo 


oo 


12) 


GT512mM10 


(si) 


512 


mix 


8.2x10 


-22 


8.0x10° 


6.2 xlO 3 


2.1 





oo 


1.1 


4.5 


11. 


0,1 


oo 


oo 


13) 


GT512mM10Bl 


(si) 


512 


mix 


8.2x10 


-22 


8.0x10° 


6.2 xlO 3 


2.1 


1 


8.2 


1.1 


5.4 


10. 


0.4 


2.8 


12. 


11) 


GT512mM10 


(s2) 


512 


mix 


8.2x10 


-22 


8.0x10° 


6.2 xlO 3 


2.2 





oo 


1.2 


8.4 


11. 


0.4 


oo 


oo 


15) 


GT512mM10Bl 


(s2) 


512 


mix 


8.2x10 


-22 


8.0x10° 


6.2xl0 3 


2.2 


1 


8.2 


1.2 


9.5 


11. 


0.4 


1.8 


10. 


16) 


GT256mM10 


(s3) 


256 


mix 


8.2x10 


-22 


8.0x10° 


6.2 xlO 3 


2.0 





oo 


1.0 


5.9 


10. 


0,1 


oo 


oo 


17) 


GT512mM10 


(s3) 


512 


mix 


8.2x10 


-22 


8.0x10° 


6.2 xlO 3 


2.0 





oo 


1.0 


5.9 


10. 


0,1 


oo 


oo 


18) 


GT512mM10Bl 


(s3) 


512 


mix 


8.2x10 


-22 


8.0x10° 


6.2 xlO 3 


2.0 


1 


8.2 


0.97 


6.4 


9.9 


0,1 


3.6 


13. 


19) 


GT256mM10B3 


(83) 


256 


mix 


8.2x10 


-22 


8.0x10° 


6.2 xlO 3 


1.8 


3 


0.92 


0.81 


8.4 


9.0 


0,1 


0.20 


2.9 


20) 


GT512mM10B3 


(S3) 


512 


mix 


8.2x10 


-22 


8.0x10° 


6.2 xlO 3 


1.8 


3 


0.92 


0.83 


8.7 


9.1 


0,1 


0.18 


2.7 


21) 


GT256mM10B10 


(S3) 


256 


mix 


8.2x10 


-22 


8.0x10° 


6.2 xlO 3 


1.8 


10 


0.08 


0.79 


6.6 


8.9 


0,1 


0.04 


1.3 


22) 


GT128cM10 




128 


comp 


8.2x10 


-22 


8.0x10° 


6.2 xlO 3 


1.8 





oo 


0.81 


1.2 


9.0 


1 


oo 


oo 


23) 


GT256cM10 




256 


comp 


8.2x10 


-22 


8.0x10° 


6.2 xlO 3 


1.8 





oo 


0.85 


1.1 


9.2 


1 


oo 


oo 


24) 


GT512cM10 




512 


comp 


8.2x10 


-22 


8.0x10° 


6.2xl0 3 


1.9 





oo 


0.87 


1.1 


9.4 


1 


oo 


oo 


25) 


GT256sM20 




256 


sol 


2.1x10 


-22 


3.2 xlO 1 


9.9 xlO 4 


4.1 





oo 


1.0 


11. 


20. 


1/3 


oo 


CO 


26) 


GT256mM20 




256 


mix 


2.1x10 


-22 


3.2X10 1 


9.9xl0 4 


4.2 





oo 


1.1 


4.5 


21. 


0,1 


oo 


oo 


27) 


GT256cM20 




256 


comp 


2.1x10 


-22 


3.2 xlO 1 


9.9 xlO 4 


4.0 





oo 


1.0 


0.60 


20. 


1 


oo 


CO 


28) 


GT256sM50 




256 


sol 


3.3x10 


-23 


2.0 xlO 2 


3.9 xlO 6 


10. 





oo 


1.1 


12. 


52. 


1/3 


oo 


CO 


29) 


GT512sM50 




512 


sol 


3.3x10 


-23 


2.0xl0 2 


3.9xl0 6 


10. 





CO 


1.1 


13. 


52. 


1/3 


oo 


oo 


30) 


GT256mM50 




256 


mix 


3.3x10 


-23 


2.0xl0 2 


3.9x10 s 


10. 





oo 


1.0 


7.0 


51. 


0,1 


oo 


oo 


31) 


GT512mM50 




512 


mix 


3.3x10 


-23 


2.0xl0 2 


3.9x10° 


10. 





oo 


1.1 


7.4 


51. 


0,1 


oo 


oo 


32) 


GT256cM50 




256 


comp 


3.3x10 


-23 


2.0xl0 2 


3.9x10° 


9.8 





oo 


0.95 


0.54 


49. 


1 


oo 


oo 


33) 


GT512cM50 




512 


comp 


3.3x10 


-23 


2.0 xlO 2 


3.9x10° 


9.9 





oo 


0.99 


0.56 


50. 


1 


oo 


oo 


34) 


GT1024cM50 




1024 


comp 


3.3x10 


-23 


2.0 xlO 2 


3.9x10° 


10. 





oo 


1.00 


0.55 


50. 


1 


oo 


CO 



Notes. Column (1): simulation name. Columns (2-10): maximum grid resolution in one direction of the three-dimensional, cu- 
bic domain, mode of forcing (solenoidal, mixed, compressive), mean density, linear box size, total mass, velocity dispersion on the box 
scale, mean magnetic-field strength (in the z-direction of the domai n), i nitial plasma /3o, and virial parameter based on Equation 1151 ■ 
Columns (11-15): time-averaged virial parameter based on Equation \16\ , computed directly from the three-dimensional gas distribution, 
the sonic Mach number, forcing parameter, ratio of thermal to magnetic pressure (plasma 0), and Alfven Mach number. To guide the eye, 
horizontal lines separate models with different sonic Mach number. 

contrast, we obtain highly inhomogeneous density distri- 
butions in our compressible, turbulent clouds. We thus 
prefer to compute a V ir based on the three-dimensional 
density field as explained abov^j In analogy, the sonic 
and Alfven Mach numbers, as well as (3 are computed as 
spatial root-mean-squared averages over all cells in the 
simulation box as a function of time, followed by aver- 
aging over time. We will show in the next section that 
this approach is justified because we find that all those 
parameters do not vary significantly with time during 
star formation. The value of the forcing parameter b was 
not determined by averaging because it was already mea- 



sured in Federrath et al. ( 2010b| Figure 8), giving best-fit 
values b = 1/3, U.4, and 1 for solenoidal, naturally-mixed, 
and compressive forcing of the turbulence, respectively. 

We do not include any data or discussion of the state 
of the clouds after SFE = 20% is reached because at that 
point in time, local feedback processes would have likely 
altered the subsequent evolution of the clouds so dras- 
tically that we cannot trust our results for higher SFE. 
Even before that, inclusion of feedback processes might 
change the results, at least locally. For example, we ex- 
pect the amount of accr eted gas to be reduced, if feed- 



back were included (e.g. 



9 Note that a similar approach is used in Herschel observations 
by | Andre et al.| l |2010| to estimate the stability of interstellar fil- 
aments, t hat is Dascd on column density instead of volume den- 
sity, but takes the spatial (projected) distribution of matter into 
account, rather than estimating the dynamical state of the cloud 
base d on the spherical, uniform-density approximation in Equa- 
tion ( ] 1 5 p . 



Wang et al.||2010| [Peters et al.| 
20111) . This fact can be accounted for by adjusting the 
local efficiency parameter e introduced in Equation Q 
to values e < 1 for all the models discussed here. We get 
back to this issue when we compare our simulations with 
the observational data in Section [6] 

The basic model parameters in Table [2 were cho- 
sen to roughly follow observed properties of molecular 
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Fcdcrrath & Klessen 



clouds, covering a range of cloud sizes L ss 0.3-200 pc, 
masses M c 300 to 4 x 10 6 , and velocity dispersions 



0.6-10 km s" 



1987 Falgarone et al 



( e.g., |Larson 1981 
'1992| 7with typica 



Solomon et al. 
cloud s earings 



summarized and discusse cT"m|M ac Low fc Klessen ( 2004 1 

However, even tnougl 



and |McKee fc Ostriker| ( |2007| 
most real clouds may roughly follow such an average scal- 
ing, the scatter around that average is typically about 
an order of magnitude or more in terms of mass, d ensity. 



and velocity dispersion for a given clou d size (e.g., Heyer 
et al.||2009| |Roman-Duval et al.||2010[ ). The procedure 
used here to determine the initial cloud parameters in 
the simulations is as follows. First, for a given target 
Mach number, we determine the appropriate size of the 
cloud by inverting the obse rve d velocity dispersion-size 
relation given by Equation (12 1. Having the size and ve- 
locity dispersi on, we then set the virial parameter given 
by Equation ( 15 ) to a value close to unity. The only 
exceptions are the M. ~ 3 models, where we set it to 
a V ir,o ~ 0.07 because this turned out to give actual virial 
parameters a v i r closer to unity after the turbulence had 
been fully established (compare columns 10 and 11 in 
Table [2]). Using the initial guess of a V i rj0 , we then solve 
for the mass of the cloud, by inverting Equation (15). 
From the mass and size, we compute the mean density 
of the model cloud. 

It is important to note that the actual virial param- 
eter obtained after two turbulent crossing times can be 
up to an order of magnitu de d ifferent from the initial 
guess provided by Equation ( 15 ), depending on the Mach 
number and forcing of the model (see Table [2 ) . This is 
because the density distribution in the state or fully de- 
veloped supersonic turbulence is highly inhomogeneous 
and is not well described by Equation ( 15 ). Thus, we do 



not know the virial parameter that arises in the regime 
of fully developed turbulence a priori. The a vlr in the 
turbulent phase is typically higher (except for the com- 
pressive forcing cases at high Mach numbers, A4 ~ 20 
and 50) than the one computed from Equation (15), 
also because we use periodic boundary conditions. This 
reduces the gravitational binding energy of the system 
com pare d to an isolated system (as assumed in Equa- 



tion 115). Real clouds are neither periodic nor isolated, 
but using periodic boundaries, we mimic the effects of the 
surrounding medium on the region studied in our compu- 
tational boxes (discussed further in Section [7]). We em- 
phasize that the virial parameters obtained here are con- 
sistent with observations, given that observational esti- 
mates of a v i r are usually obtained based on Equation ( 15 1 
or column- density versions of it. 

Magnetic-field strengths for the MHD simulations were 
chose n to be consistent with the range obse rved in clouds 
(e.g., |Crutcher|1999[[Crutcher et al.||2010[ ). We vary the 
magnetic held tor simulations with mixed forcing and 
fixed sonic Mach number of A4 ~ 10, which gives us a 
good indication of the role of magnetic fields for ty pical 
molecular cloud prope rties. Heiles & Troland (2005) and 
Crutcher et al. (2010) show that most clouds with num- 
ber densities in the range 10-10 4 cm -3 have magnetic- 
field strengths in the range B z w 1-10 /iG, with an ap- 
parent peak of the distribution at around B z « 3/iG. 
Our MHD simulations have mean densities of about 
200 cm~ 3 , so we decided to compare models with B z — 1, 



3, and 10 /xG, in order to cover the observed range of line- 
of-sight magnetic-field strengths. 

4. SIMULATION RESULTS 

After the initial turbulent state has been e stab lished by 
driving for two crossing times (see Section 3.4) in each 
simulation, we study the subsequent evolution under the 
influence of self-gravity by looking at column density pro- 
jections of the simul ated clouds and their magnetic-field 
morphology (Section |4.1[ ). We then discuss the time evo- 
lution of Ovir, M., M.a, and SFE and measure the SFR 
in Section |j~2l 



4.1. Cloud and Magnetic- field Morphology 
4.1.1. Effects of the Magnetic Field 

Figure [4] shows the time evolution of column density 
snapshots (from top to bottom) for models with mixed 
forcing at M. = 10 and 512 3 resolution for initial mag- 
netic fields i?o = 0, 1, and 3^G (left, middle, and right 
panels). Key initial parameters (box size, total mass, 
etc.), the time in units of tg(po), the SFE, and the num- 
ber of sink particles formed are given in each panel. The 
top row shows the gas at t — 0, i.e., when turbulence is 
fully developed and self-gravity is switched on. We see 
shocks and large-scale structure induced by the large- 
scale turbulence with column density contrasts ranging 
over more than four orders of magnitude. Comparing 
the purely HD run (left) with the two magnetized runs 
(middle and right), we see that shocks become smoother 
and density contrasts slightly decrease as the magnetic- 
field strength increases. This is because magnetic fields 
act like a cushion, reducing density fluctuations, due to 
the additional magnetic pressure parameterized either by 
plasma f3 or the Alfven Mach num ber A^a (see Equa- 
tion [|or[5| and |Molina et al.||2012[ ), the time -averaged 
values of which are given in Table [2j At later times, 
the gas starts collapsing locally at sites previously com- 
pressed by the supersonic turbulence, at which point lo- 
cal filaments become more and more massive as they ac- 
crete gas from the surrounding and eventually become 
so dense that these cores have to be replaced with sink 
particles, allowin g us to advance the simulations to later 
times (see Section 3.3 ). The radius r s j n k of the sink parti- 
cles is determined by the numerical resolution constraint, 
and is given in each panel, as soon as sink particles have 
formed. Our resolution is insufficient to resolve individ- 
ual stars, but the sink particles can be regarded as dense, 
bound cores in our simulations. 

Comparing the runs with different magnetic-field 
strengths in Figure HI we see two important effects with 
increasing magnetic field: (1) a reduction of fragmenta- 
tion, i.e., fewer sink particles have formed by the end of 
the simulations at SFE = 20% and (2) reaching a given 
SFE takes longer, i.e., the core formation rate and hence 
the SFR are reduced. For instance, when the SFE has 
reached 20%, the runs with Bq = 0, 1, and 3/iG have 
formed 109, 71, and 63 sink particles in 0.71, 0.85, and 
1.1 tff (po) , respectively. 

The higher the magnetic field, the larger the 
topologically-connected structures, compared to the 
more fragmented and dispersed filaments in the purely 
hydrodynamical run. Comparing numerical simula- 
tions and observations with filament-tracking tools (e.g., 
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Mach 1 




log i0 ( column density [g/cm ] ) 



2.8 -2.6 -2.4 -2.2 -2.0 
log, ( M s]nk / M 

box / 



Fig. 4. — Time evolution of column density projections of the simulations with mixed forcing at Ai = 10 with initial magnetic- field 
strengths Bo = /J.G (left panels), Bq = 1 fj,G (middle panels), and Bo = 3 fiG (right panels). The times shown correspond to the initial, 
fully developed turbulent state, t = (top panels), and when the star formation efficiency reached SFE = 5% and 20% (middle and bottom 
panels, respectively), i.e., 5% and 20% of the gas was accreted by sink particles (shown as circles with the sink particle radius). The higher 
the magnetic field, the slower the star formation (see time in the top right corner), and the fewer sink particles form (see bottom right 
corner of each image) due to the increasing magnetic pressure. 
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Fedcrrath & Klessen 













0,0 





Mach 10, mix, 5 = 3^G 
z- projection 



= 0.75 




Fig. 5. — Column density and magnetic-field vectors in simula- 
tion model #20 (GT512mM10B3; see Tablepl with B = 3/iG for 
t/t ff = (top) and SFE = 10% (bottom). T'he magnetic field is 
amplified in the core and cluster regions, where compression and 
turbulent dyna mo action both contribute to increasin g the field 
strength locally | |Sur et al,|2 010 Fedcrrath ct al. 2011c|. The mag- 
netic field frequently changes direction and strength m the cores 
in this model of super-Alfvenic turbulence with Ma ~ 2.7 (see 
Table pi, best seen in the movies (see additional online material) . 
The coTors are identical to Figure H] 



Andre et al.|[20T0j |Men'shchikov eFaLl[20T0l |Arzouma- 
nian et al.||2011| IHill et al.||2011| ISchneider et al.||2012| 



or polarization analyses (e.g., Burkh art et al.||2(J12 ) may 
eventually help to reveal the role of magnetic fields. In 
particular, the orientation of magnetic fields might tell us 
about its dynamical influence (ISchneider et al.]2010" Li fe 
Henning|20lT| |Peretto et al.|2012| . In Figure|5} we snow 
the time evolution of column density snapshots with lo- 
cal magnetic-field vectors computed by a mass-weighted 
average along the line of sight superimposed, for the run 
with B = 3/iG for t/t s = (top) and SFE = 10% (bot- 
tom) . The magnetic field grows due to compression of 
the field lines and due to dynamo action (Sur et al.|2010[ 
Federrath et al.|2011c Bertram et al.|20T2 ), particularly 
in regions where dense cores accumulate and form clus- 
ters. The magnetic field is very intermittent and shows 
no particularly preferred direction in the cluster centers 
because the gas motions are so chaotic that the magnetic- 



field direction changes frequently. The magnetic field is 
of moderate strength compared to the turbulence in this 
case, shown by the average super-Alfvenic Mach num- 
ber in this simulation, Ma ~ 2.7 (see Table [2]). The 
field strengths are consistent with observations in typical 
molecular clouds. On scales larger than molecular clouds 
and on Galactic scales though, the turbulence might be 
trans- Alfvenic rather than super-Alfvenic, which would 
naturally le ad to more aligned magnetic field structu res 



there (e.g.,|Beck et al 
fe Henning|[201l| ) 



1996 Heiles & Troland 2005 Li 



4.1.2. Effects of Turbulent Forcing and Sonic Mach Number 

After having looked at the time evolution of column 
density snapshots in mixed, Ai ~ 10 simulations, we 
now focus on the gas morphology when SFE = 10%, rep- 
resenting a typical molecular cloud value, but compar- 
ing different forcing and sonic Mach numbers. Figure [6] 
shows column density projections of the 512 3 runs with 
solenoidal forcing (left panels) and compressive forcing 
(right panels) at M. ~ 3 (top), M. ~ 10 (middle), and 
Ai ~ 50 (bottom). Note the different length and mass 
scales probed in these images, with box sizes of L = 0.3, 
8 and 200 pc, and masses of M c = 310, 6.2 x 10 3 and 
3.9 x 10 6 M Q , respectively. Since the resolution is fixed, 
the sink particle radii vary from r S j n k = 335 AU over 
fsink = 0.04 pc, up to 1 pc. Thus, neither of those rep- 
resents stars, but rather star clusters in the largest-scale 
runs and potentially protostellar accretion envelopes in 
the smallest-scale runs. The scale and mass sequence 
from the bottom to the top panels in Figure [6] can be in- 
terpreted as zooms into patches of larger-scale runs and 
re-simulating these patches with higher resolution in suc- 
cessively smaller boxes. Clearly, these images emphasize 
how artificial this kind of numerical experiment is, yet 
real molecular clouds exhibit similar hierarchical struc- 



tures ( |Falgarone et al.|1992||Osse nkopf & Mac Low|2002 
often characterized as fractals ([Scalo 1990 Elmcgrccn 



Falgaronc 1996; S tutzki et all p 
Roman-Dtival et al.||2010p . Th 



1 



0981 iSanctiez'eFar 2005 ; 



e fractal dimension D in- 
terred from different techniques (A-variance, box count- 
ing, mass-size relation, and perimeter-area method) was 
shown to vary between D « 2.6 for purely solenoidal and 
D ss 2.3 for purely compressi ye forcing, in the ran ge of 
observational determinations ( Federrath et al. 2009 ) , and 
is consistent with theoretical ideas to explain t he slope 
of the stellar IMF (IChabrier fc Hennebelle|2011[ ). As can 
be seen in Figure |oJ compressive forcing produces more 
sheet-like structures (planar shocks), while solenoidal 
forcing produces more volume-filling structures, provid- 
ing a visual explanation for the dependence of D on the 
forcing. 

Besides the morphological distinctions, the most strik- 
ing difference between solenoidal and compressive forc- 
ings is the timescale of core and star formation (compare 
i/tff in the upper right corner of each panel in Figure [6]). 
For fixed Mach number, cloud size, and mass, compres- 
sive forcing accelerates the conversion of gas into stars 
compared to solenoidal forcing by factors of 4, 8, and 
12 for the M. ~ 3, 10, and 50 runs, respectively, when 
SFE ~ 10%. This result emphasizes the important role 
of the turbulent forcing for setting the SFR. 

4.2. Time Evolution of a v [ v , Ai, A4a, and SFE 
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Fig. 6. — Column density projections of the simulations with solenoidal forcing (left panels) and compressive forcing (right panels) for 
Mach numbers Ai ~ 3 (top), Ai ~ 10 (middle), and ~ 50 (bottom), when 10% of the initial gas mass is accreted by sink particles 
(shown as circles with the sink particle radius). The mass and size of the three-dimensional domains, and the number of sink particles 
formed, are given in each panel. In addition to the morphological differences between the forcings for a given Mach number, the elapsed 
time in units of the frccfall time at the mean density (see label in the top right corner of each panel) is significantly different between the 
two extreme cases of turbulent forcing, suggesting extremely different star formation rates for solenoidal and compressive forcing. 
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Fig. 7. — Time evolution of the rms Mach number A4 (top), the 
virial parameter a VIr (middle), and the star formation efficiency 
SFE (bottom) for models with fixed Mach number Ai ~ 10, and 
without magnetic field, but with solenoidal forcing (gold) and com- 
pressive forcing (violet) at numerical resolutions of 128 3 (dotted), 
256 3 (dashed), and 512 3 (solid), as well as for three mixed-forcing 
models (black), each forced with a different random number se- 
quence at fixed 512 3 resolution: seedl (dash-dotted), seed2 (triplc- 
dot-dashed), and seed3 (solid). The time axis is scaled in units of 
the freefall time at the mean density of the respective simulation 
(see Table . Values of SFRg , measured from linear fits to the 
SFE-time curves (bottom panel) in the range SFE = [4%, 20%] 
are given for each model in the legend. 

4.2.1. Effects of the Forcing, Random Seed, and Resolution 

We now turn to the detailed time evolution and deter- 
mination of the SFR in the simulations. Figure [7] shows 
the time evolution of dynamical quantities M., cx v i T , and 
SFE for models with M. ~ 10 and Bo = for solenoidal 
and compressive forcing at numerical resolutions of 128 3 , 
256 3 , and 512 3 grid cells, and for mixed forcing with three 
different random seeds of the turbulent forcing. The 
Mach number (top panel) shows variations of order 10% 
around the target Mach number of A4 ~ 10 for each 



simulation, and some systematic variations with differ- 
ent forcing and random seeds. The differences between 
solenoidal, mixed, and compressive forcings are caused 
by stronger dissipation with more compressive forcing, 
requiring a higher forcing amplitude to reach the same 
Mach number than in solenoidal forcing. We adjust the 
amplitude of the forcing such that the gas reaches a given 
Mach number in the fully developed turbulent phase. 
Since the value of M. depends on nonlinear dissipation 
properties, i.e., strengths of shocks and amount of vor- 
ticity gener ated, and thus on the M ach number of the 
turbulence ( |Federrath et al.||2011b| ), the time-averaged 
rms Mach number tor a given forcing amplitude cannot 
be predicted a priori and must be adjusted iteratively by 
running test simulations with different forcing amplitude 
and measuring the time-averaged rms Mach number, re- 
sulting in some deviation of the actual Mach number 
from the target Mach number (see the time-averaged M. 
for each model in Table [2]) . The temporal fluctuations 
and the differences between random seeds, however, are 
purely statistical. In order to compare our simulation 
data with the analytic theories, we thus always use the 
volume- and time-averaged quantities entering the theo- 
retical models from Section [2.51 

The middle panel of Figure [7] shows a v -„(t). As for 
M(t), the resolution dependence is only marginal, and 
significantly less than t he statistical fluctuations (see also 
Kitsionas etaLl[2fJ09l |Price fc Federrath|[2010l |Kritsuk| 



et al.||2011b| showing that one-point statistics are typi- 
cally well converged with grid resolutions of 256 3 cells). 
This demonstrates that the length scales of the dominant 
gravitational structures are resolved well enough in our 
present numerical experiments and that our definition of 
a v ir is robust with respect to changes in the numerical 
resolution. 

The difference of a v ; r = 2 -Ekin/l-E-gravl between the 
forcings deserves some attention. While the A4 ~ 10 
runs with solenoidal forcing have a v i r ~ 12, the com- 
pressive ones have a v h ~ 1.1 (see Table [2]), even though 
the Mach number is similar and the mass of the clouds 
is identical. In Figure [6] we saw that compressive forc- 
ing produces more locally compressed structures than 
solenoidal forcing, resulting in an overall higher gravi- 
tational binding energy |-E gra v| compared to solenoidal 
forcing. The total kinetic energy E'kin on the other hand 
is the same within a factor of ~ 2, which means that 
the factor of ~ 10 difference in a V i r is primarily due 
to the difference in |-E grav |. This shows that compar- 
ing simple theoretical estimates of the virial parameter, 
solely based on the total mass as a measure for |-E grav | 



(as, e.g., assumed in Equation [15]) , should be considered 
with great caution because such an estimate ignores the 
internal structure of the clouds. Thus, we prefer to es- 
timate a v i r based on the actual spatial distribution, as 
we have done in Figure [7] and in Table [2] for all models. 
We emphasize that this direct comparison of a V i r ,o with 
a v ir performed here means that observational estimates 
of the virial parameter b ased on global measures such 
as described by Equation ( 15 ) or alike are only accurate 
within an order of magnitude. Measurements of gravi- 
tational (in) stability based on the actual column density 
distribution of fila ments in Herschel o bservations of the 
Gould Belt (GB) by |Andre et al.| ( |2010"| ), for example, are 
thus likely more accurate and meaningful than estimates 
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based on uni form -density, spherical approximations such 
as Equation ( 15 1. 

The bottom panel of Figure [7] shows the time evolution 
of the total mass accreted by sink particles, divided by 
the total cloud mass, i.e., the SFE. We measure the slope 
of these curves by fitting a linear function in the interval 
SFE = [4% , 20%], which gives the SFR ff for each model 
quoted in the legend. We choose to set the lower limit of 
the fit range to SFE = 4% because the initial accretion 
phase is highly nonlinear with a fast increase in slope, 
after which the accretion becomes roughly linear in time, 
such that the slope is reasonably well defined for most of 
the models. 

First, we study the dependence of SFRff on the random 
seed. The three models with mixed forcing and different 
random seeds (seedl, seed2, seed3) exhibit variations in 
SFRff by a factor of 1.5. However, other seeds might de- 
viate further from this, such that the factor 1.5 in SFRff 
is a lower limit for the uncertainty introduced by the 
random seed. When we compare mixed-forcing models 
at M ~ 10 with different magnetic-field strengths later, 
we always compare runs with seed3 because the SFRff 
for seed3 is in between the ones measured for seedl and 
seed2, thus giving the best average behavior for the data 
at hand. 

Finally, we investigate the resolution dependence of 
our simulations with solenoidal and compressive forc- 
ings in Figure For resolutions of 128 3 , 256 3 , and 512 3 
grid cells, we find that SFR ff = 0.27, 0.17, and 0.14 for 
solenoidal forcing, and SFR ff = 1.58, 2.27, and 2.75 for 
compressive forcing, respectively. Thus, with increas- 
ing resolution, SFRff is decreasing for solenoidal forcing, 
but increasing for compressive forcing. The difference 
in SFRff between 128 3 and 256 3 is a factor of 0.63 for 
solenoidal forcing, and a factor of 1.44 for compressive 
forcing. These factors become smaller when we com- 
pare the 256 3 with the 512 3 simulations, giving factors 
of 0.82 for solenoidal forcing and 1.21 for compressive 
forcing. Thus our results converge with increasing res- 
olution. Moreover, we can estimate SFRff in the limit 
of infinite resolution from extrapolating the convergence 
behavior. Doing this, we see that our measurements of 
SFRff at 128 3 resolution are converged only within a fac- 
tor of about 2.5, so we discard the two 128 3 simulations 
(GT128sM10 and GT128cM10) in all the following. In 
contrast, the 256 3 data are converged to within a fac- 
tor of 1.5 for both solenoidal and compressive forcings, 
which is similar to the uncertainty introduced by varying 
the random seed as discussed in the previous paragraph. 
Thus, differences larger than a factor of 1.5 in SFRff 
between models with different physical parameters are 
likely of physical rather than numerical or statistical ori- 
gin. For instance, the SFRff for compressive forcing with 
A4 10 is more than an order of magnitude larger than 
the SFRff of the respective solenoidal simulation, demon- 
strating the physical importance of the turbulent forcing 
for controlling the SFR. 

4.2.2. Effects of Increasing the Sonic Mach Number 

Having looked at models with different forcing, random 
seed, and resolution, we now study models with varying 
Mach number. Figure [8] shows the same as Figure [7j but 
for compressive-forcing models with Mach Ai ~ 3, 5, 10, 
20, and 50. The Mach number increases slightly with 
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Fig. 8. — Same as Figure [7] but for compressive-forcing models 
with sonic Mach number, A4~ 3, 5, 10, 20, and 50. 

time in all models, which is caused by local accelerations 
during collapse. This only accounts for a few percent at 
most. The exception is the A4 ~ 3 model, for which M. 
increases by almost a factor of two (top panel). This is 
because the M. ~ 3 model is more gravitationally unsta- 
ble initially, indicated by the virial parameter (middle 
panel), which increases to around unity, similar to the 
other models. The SFRff generally increases with Mach 
number due to the stronger local compressions created at 
higher M. (bottom panel). The only exception is again 
the M. ~ 3 model, which has a slightly higher SFRff 
than the A4 ~ 5 model because the M. ~ 3 model is 
more unstable and starts collapsing globally, while this 
is not the case in the other models. The difference of 
SFRff between M. ~ 5 and 50 is about a factor of 4.4. 

4.2.3. Effects of Increasing the Magnetic- field Strength 

Finally, in Figure [9] we investigate the time evolution 
of models with different initial magnetic-field strengths, 
B = 0, 1, 3, and 10 ^G (initial plasma /3 = 8.2, 0.92, 
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Fig. 9. — Same as Figure [7] but for mixed-forcing models (6 = 
0.4) at M ~ 10 with diffcrentinitial magnetic-field strengths Bq = 
0, 1, 3, and 10 fiG. The additional panel on the top shows the time 
evolution of the Alfven Mach number in the MHD simulations. 

and 0.082; see Table [2]). The panels are the same as in 
FigurelTJ except for an additional panel on the top, show- 
ing the^lfven Mach number Ma- Apart from some tem- 
poral fluctuations, Ma, M, and a v ; r are fairly constant 
over time. Both Ma and M show some minor systematic 



decrease, which is caused by dynamo action, amplifying 
the magnetic fiel d by converting turbulent energy int o 
magnetic energy (Brandenburg & Subramanian 2005). 
Most of the dynamo action, however, took place already 
during the first two turbulent crossing times, t < Otg, 
during which the turbulence becomes fully established 
(compare columns 9 and 14 of Table |2| . The dynamo is 
nearly saturated at t — with only very slow linear am- 
plification happening afterward. In addition, field lines 
are compressed during local collapse, amplifyingthc field 
further in dense cores and clusters (see Figure pi. 

Most importantly, the last panel of Figure [9] snows that 
the SFRff decreases monotonically with increasing mag- 
netic field because of the stabilizing effect of the magnetic 
pressure. The strongest magnetic field case studied here 
(B = 10 fiG, M A ~ 1.3) has an SFR ff w 0.24, which 
is almost a factor of two smaller than in the respective 
purely hydrodynamical run (B = 0, SFRff ss 0.46). A 
similar reduction of the SFR with strong magnetic fields 
compared to purely hydro dynamical or weakly magne- 
tize d models is reported in Padoan & Nordlund (20111 



and Padoan et al. 
by a factor ot 



(2012), who find a maximum reduction 
o! This is a significant, but relatively 
small effect compared to the influence of different forc- 
ing on the SFRff (see above). Magnetic fields reduce 
SFRff, but are unlikely the major player in controlling 
the SFR, provided that molecular cloud turbulence is 
super-Alfvenic or at most trans-Alfvenic. This seems to 
be the case in most clouds. However, as pointed out ear- 
lier, on larger scales than molecular clouds, i.e., in the 
warmer, mainly atomic part of the IS M, turbulence may 
be trans-Alfvenic or even sub-Alfvenic B eiles fc Trcjand 
20051 |Li & Henningl 120111 iHeyer & Brunt||2012l), ren 



dering magnetic fields potentially more important in the 
process of molecular cloud formation. Still, even inside 
molecular clouds, magnetic fields seem to reduce frag- 
mentation significantly (see Figure EJ), thus potentially 
having a strong imp act on the mass distr i bution of cores 
and 
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5. COMPARING SFRS IN THE MHD SIMULATIONS WITH 
THEORETICAL PREDICTIONS 

Using the dimensionless parameters a V i r , M, b, and 
P (or A^a) measured for each numerical simulation and 
listed in the last five columns of Table [2] we can now 
compute the SFRff predicted by each of the six theories: 
KM, PN, HC, and multi-freefall KM PN, HC, introduced 
in Section [2] (summarized in Table [Tj) , and compare it to 
the simulated SFRff. The comparison between SFRff 
(theory) and SFRff (simulation) is shown in Figure 10 
(left panels: KM, PN, HC; right panels: multi-freefalT 
KM, PN, HC). The SFR ff in each of the six theoretical 
models is fully determined by a V i r , M, 6, and /?, except 
for the parameters e/(j) t and the fudge factors 4>x (KM), 
9 (PN), or y cut (HC). In the simulations, the local ef- 
ficiency e = 1 because we did not include any form of 
feedback, but l/4>t and the theory fudge factors are free 
parameters. In order to constrain them for each the- 
ory, we perform two-parameter fits of SFRff (theory) to 
SFRff (simulation). The best-fit parameters are listed 
in the legend of Figure 10 Table [3] additionally lists 



uncertainty estimates for the parameters, together with 
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Fig. 10. — SFRff (theory) for the six theories listed in Table [T] KM (boxes), PN (diamonds), and HC (crosses) in the left panels, and 
the corresponding multi-freefall versions of the theories in the right panels, computed based on the numerical simulation parameters a v ir, 
A4, b, and j3 listed in Table p] and compared with the SFRff (simulation). The simulation number is given in each of the KM boxes. The 
analytic model predictions, SFRff (theory), were fitted to SFRff (simulation) with the fit parameters e/0t (where e = 1 by definition in the 
simulations) and the fudge factors <f> x (KM), 8 (PN), and y cu t (HC). The best-fit parameters are given in the legend. The fits in the top 
panels only used the hydrodynamic models for which Bo = 0, while the fits in the bottom panels include all MHD models listed in Table [2] 
(except for the low-resolution 128 3 -models). A zoom of the region containing the MHD models is shown in the inset plots in the bottom 
panels, where only the six MHD simulations are included. The diagonal solid line in each plot represents perfect agreement between SFRff 
(theory) and SFRff (simulation). The best-fit parameters with uncertainties and x 2 -values are listed in Table [3] Each simulation-theory 
data pair is listed in Table [4] 



X 2 -values, the number of degrees of freedom (DOF) in 
the fits, and the reduced Xrcd = X 2 

/DOF. The x 2 cd is a 
quantitative indicator for the goodness of fit, with better 
fits having smaller Xrcd- To separate the effects of the 
magnetic field, we only use purely HD models (Bo = 0) 
in the top panels of Figure 10 (HD fit), while we in- 
clude all MHD models in the Dottom panels (MHD fit). 
This distinction is also made in Table [3| Inset plots in 
the bottom panels show a zoom-in on the MHD models 
only. The solid diagonal line in each panel represents 
SFRff (theory) = SFRff (simulation), i.e., perfect agree- 
ment bet wee n theory and simulation. 

Figure [TO] shows that all the theoretical models ex- 
hibit some positive correlation between SFRff (theory) 
and SFR ff (simulation). The multi-freefall KM and PN 
models (right panels) show much better agreement with 
the simulation data in both the HD and MHD fits, indi- 
cated by the smallest x 2 0C i = 1-2-1.3 (see Tablepj), than 
the original KM and PN models (left panels). The HC 
models exhibit the opposite behavior, i.e., the HC the- 
ory gives slightly better fits than the multi-freefall HC 
theory. This is not surprising because both HC mod- 
els use the multi-freefall factor, but the HC model ad- 



ditionally includes turbulent sup port in the estimate of 
the threshold density (Equations 38 and 39 into Equa- 



tion 



37), while the multi-freefall HC modcrahly includes 



thermal support (Equation 38 only). However, all HC 



fits exhibit relatively large xf cd ~ 4.9-6.2. The reason 
for this is the choice of the critical density in the HC 
models and its resulting dependence on the sonic Mach 
number, p cr ; t oc A4~ 2 , while all KM and PN models have 
Pent kM +2 . which is (apart from the different choice of 
fudge factors) the only fundamental difference between 
the multi-freefall HC and the two multi-freefall KM and 
PN models (see Table [l]). The difference in fudge factors 
is irrelevant in this comparison because they all enter 
in the same way for each theory, simply as factors in 
the critical density, for which the fitting procedure de- 
termines the best-fit value automatically. In contrast, 
the dependence of SFRff (theory) on a vlr , M., b, and /3 is 
determined by each analytic theory separately. Table [JJ 
gives an overview of the basic similarities and differences 
between the six theoretical models for the SFRff. 

The KM fits also exhibit fairly large x 2 od = 5.3 and 5.7 
in the HD and MHD fit set, respectively. In contrast, the 
multi-freefall version of the KM model gives much better 
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TABLE 3 

SFRff (Theory) -SFRff (Simulation) Fit Parameters (Fig. [Toll. 



(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


Theory {HD fit) 


V<k 


Fudge Factor 


x 2 


DOF 




KM 


3.00 ±n/a 


<f> x = 0.12 ± n/a 


127 


24 


5.3 


PN 


1.50 ± 0.16 


6 = 0.65 ± 0.05 


46 


24 


1.9 


HC 


0.24 ± n/a 


Vent = 1-3 ± n/a 


135 


24 


5.6 


multi-ff KM 


0.49 ±0.06 


4> x =0.19 ±0.02 


32 


24 


1.3 


multi-flf PN 


0.49 ±0.06 


e = 0.97 ± 0.10 


32 


24 


1.3 


multi-ff HC 


0.21 ± n/a 


2/cut = 1-1 ± n/a 


119 


24 


6.2 


Theory (MHD fit) 


KM 


4.10 ± n/a 


<j} x = 0.17 ± n/a 


172 


30 


5.7 


PN 


1.40 ± 0.14 


6 = 0.70 ± 0.04 


54 


30 


1.8 


HC 


0.21 ± n/a 


y cut = 4.5 ± n/a 


147 


30 


1.9 


multi-ff KM 


0.46 ±0.06 


4> x = 0.17 ± 0.02 


39 


30 


1.3 


multi-ff PN 


0.47 ± 0.06 


e = i.o ±o.i 


37 


30 


1.2 


multi-ff HC 


0.20 ± n/a 


2/cut = 5.9 ± n/a 


152 


30 


5.1 



Notes. Column 1: Theoretical model according to Ta- 
ble [TJ Columns 2 and 3: Fit parameters for the HD fit set 
(top) and MHD fit set ( bot tom), corresponding to the top and 
bottom panels in Figure |l0| Column 4: \ 
5: Number of degrees of freedom (DOF) 
numerical models used for fitting (see Table FIl minus 2 (the 
number of fit parameters). The last column (6) snows the reduced 



2 of the fit. Column 
i.e^the number of 



X /DOF, enabling a direct comparison of the fit quality 



between the HD and MHD fit sets. Smaller X 2 ec j indicate better 
fits. Uncertainty estimates for the fit parameters in columns 2 and 
3 are only shown for models with X 2 e d <2. 

fits (Xr C d = 1-3 for both the HD and MHD fits, respec- 
tively). The original PN model already gives fairly good 
fits (x 2 rcd = 1-9 and 1.8), but again, the multi-frccfall PN 
version gives better fits, in fact the best fits of all ana- 
lytic theories (x? e d = 1-3 f° r the HD and Xrod = 1-2 for 
the MHD fit). The HD fits for the multi-freefall KM and 
multi-freefall PN models are identical because in the HD 
limit the two theories are identical, while in the MHD 
case, the only difference is the /^-dependence of p r , it , 
which is p cr it,KM oc 1/(1 + P^ 1 ) for KM (Equation 20), 
while it is p C rit,PN oc /(/3) given by Equation (31) for 



the PN theory. However, the difference in x^ed between 
multi-ff KM and multi-ff PN is very small, such that both 
the multi-freefall KM and multi-freefall PN models pro- 
vide the best match to our set of numerical simulations. 

The best-fit MHD theory parameters for the multi- 
freefall KM and multi-freefall PN models are similar (see 
Table [3]). Taking into account the full range of error 
margins, we find l/</> t = 0.4-0.55, and 4> x — 0.15-0.21 
and = 0.87-1.1. The multi-ff KM fit thus suggests 
a close correspondence of the magnetothermal Jeans 
length (Equation [2l|) and the magnetosonic scale (Equa- 
t ion [22]) with a correction of order <f> x = 0.18±0.03. The 
multi-It PN model fit supports the expected large-scale 
injection of turbulence, parameterized by 9 = 0.99 ±0.11 
(see Section [2]). Moreover, the \rcd = 1-2-1-3 of the 
multi-ff KM and multi-ff PN fits are similar, but slightly 
smaller in the MHD fit set than in the HD fit set. This 
indicates that the magnetic-field dependence in the an- 
alytic models provides a good match to the simulation 
data, and that our ext ension of the multi-ff KM model 
to MHD in Section l2~4l is reasonable. 

Even though the agreement between SFRff (theory) 
and SFRff (simulation) is very good for the multi-ff KM 
and multi-ff PN models shown in Figure \10\ some nu- 



merical simulations only agree within a factor of 2-3 
with the analytic prediction. To distinguish each sim- 
ulation, we added the sim ula tion numbers of Table [2] 
in each KM box of Figure [lO] The values of the mea- 
sured SFRff (simulation) ana the computed SFRff (the- 
ory) are listed in Table El Generally, the multi-ff KM 
and PN theories agree with the simulation data within a 
factor of two. The simulation with the largest deviation 
is model #30 (GT256mM50), for which the predicted 
SFRff by the multi-ff KM and PN models is a factor of 
2.9 and 2.7 higher than the measured SFRff in the sim- 
ulation. The higher-resolution version of this simulation 
with 512 3 cells (#31: GT512mM50) shows an improve- 
ment, such that SFRff (simulation) is now only a factor 
of 2.2 higher than SFRff in both the multi-ff KM and 
PN theories. A similar trend with increasing resolution 
is obtained for MHD models #19 (GT256mM10B3) and 
#20 (GT512mM10B3), as well as for #28 (GT256sM50) 
and #29 (GT512sM50), all converging toward the diago- 
nal, solid line in Figure [10] for the multi-freefall KM and 
PN models. This improvement with increasing resolution 
can be seen best for the M. ~ 50, compressive-forcing 
models #32 (GT256cM50) with 256 3 , #33 (GT512cM50) 
with 512 3 , and #34 (GT1024cM50) with 1024 3 resolution 
in the right panels of Figure |10| The convergence with 
increasing resolution suggests that the analytic theories 
give reasonable results and that we have constrained the 
theory parameters well with our set of numerical simula- 
tions. 

The overall agreement between the theories and simu- 
lations is encouraging. Although some numerical models 
only agree within a factor of 2-3 at the limited resolu- 
tion available, we have to keep in mind that the overall 
agreement holds over two orders of magnitude in SFRs, 
from SFRff ss 0.1 to 10, as covered by all the numerical 
simulations with different virial parameters, Mach num- 
bers, f orci ng, and magnetic-field strengths, combined in 
Figure |T0] All our simulations are fit simultaneously by 
the muitl^ff KM and multi-ff PN models. 

6. COMPARISON WITH OBSERVATIONS 

Here we compare the MHD simulation results of the 
SFR from Section [5] with observations of Galactic clouds. 
Since observed SFRs are usually quoted as SFR column 
densities, EsfRj i-e., SFR per unit area, we convert the 
simulated SFRs to Ssfr to facilitate the comparison with 
observations. 

6.1. MHD Simulations Converted to S gas and Ssfr 

We measure E gas and Esfr with a method as close 
as possible to wha t observers do to infer SsFR.-to-S 
relati ons (see, e.g., |Bigiel et al 



2008 Heiderman et al. 



2010 1, including the ettects of telescope beam smooth- 



mg. For each simulation, we construct two-dimensional 
projections of the gas column density E gas and the sink 
particle column density Sgp along each coordinate axis: 
x, y, z. All maps were smoothed to a resolution N ICS /8 
with the numerical resolution N res given in Table [2j 
such that the size of each pixel in the smoothed maps 
slightly exceeds the si nk p article diameter (which is 5 
grid cells; see Section 3.3). We also test smoothing to 
iVre S /32 below, which yields similar results. We then 
search for pixels with a sink particle column density 
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TABLE 4 

SFRff in the Simulations Listed in Tableland Theoretical Predictions for the Best-fit MHD Parameters in Table [3] 



Model 


SFR ff : 


Simulation 


KM 


PN 


HC 


multi-ff KM 


multi-ff PN 


multi-ff HC 


(1) 






(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


01) 


GT256sM3 




6.2X10" 1 


3.4x10+° 


7.6X10" 1 


2.7X10" 1 


5.3X10" 1 


5.3X10" 1 


2.6X10" 1 


02) 


GT512sM3 




6.2X10" 1 


3.3x10+° 


7.4xl0 _1 


2.7X10" 1 


5.3X10" 1 


5.3X10" 1 


2.6xl0 _1 


03) 


GT256mM3 




7.3X10" 1 


3.5x10+° 


9.1X10- 1 


3.0X10" 1 


6.1xl0 _1 


6.1X10- 1 


2.8X10- 1 


04) 


GT256cM3 




2.5x10+° 


4.0x10+° 


8.9X10- 1 


4.9X10- 1 


1.1x10+° 


1.1x10+° 


4.6xl0^ 1 


05) 


GT512cM3 




2.4x10+° 


4.0x10+° 


9.1X10" 1 


4.9X10" 1 


1.1x10+° 


1.1x10+° 


4.6X10" 1 


06) 


GT256sM5 




2.4X10- 1 


2.8X10- 1 


8.2xl0~ 2 


3.0X10- 1 


1.3X10- 1 


1.2X10" 1 


3.3X10- 1 


07) 


GT256mM5 




2.5X10" 1 


7.2X10- 1 


2.9X10" 1 


3.6X10- 1 


3.0X10- 1 


2.9X10" 1 


3.6X10" 1 


08) 


GT256cM5 




2.1x10+° 


3.0x10+° 


1.5x10+° 


6.7X10- 1 


1.4x10+° 


1.4x10+° 


6.3X10" 1 


09) 


GT128sM10 




2.7X10" 1 


n/a 


n/a 


n/a 


n/a 


n/a 


n/a 


10) 


GT256sM10 




1.7X10" 1 


1.3X10- 1 


1.4X10" 1 


4.9X10" 1 


1.6X10" 1 


1.5xl0 _1 


5.2X10" 1 


11) 


GT512sM10 




1.4X10" 1 


1.3X10" 1 


1.3X10" 1 


4.9X10" 1 


1.6xl0 _1 


1.5xl0 _1 


5.2X10" 1 


12) 


GT512mM10 


(scedl) 


5.8X10- 1 


5.9X10- 1 


6.3X10- 1 


6.2X10- 1 


5.6X10" 1 


5.5X10" 1 


6.0X10- 1 


13) 


GT512mM10Bl 


(sccdl) 


4.6X10" 1 


5.4X10" 1 


5.5X10- 1 


5.4X10" 1 


4.4X10" 1 


4.8XHT 1 


5.3X10- 1 


11) 


GT512mM10 


(seed2) 


3.9X10" 1 


3.1X10- 1 


4.0X10" 1 


6.2X10" 1 


3.9X10" 1 


3.7X10" 1 


6.2X10" 1 


15) 


GT512mM10Bl 


(seed2) 


2.9X10" 1 


2.9X10- 1 


3.4X10- 1 


5.1X10- 1 


2.8X10- 1 


3.2X10" 1 


5.2X10- 1 


16) 


GT256mM10 


(sccd3) 


4.6X10" 1 


4.6X10" 1 


4.9X10" 1 


5.9X10- 1 


4.6xl0 _1 


4.4xl0~ x 


5.8X10- 1 


17) 


GT512mM10 


(seed3) 


4.6X10" 1 


4.6X10" 1 


4.9X10- 1 


5.9X10- 1 


4.6X10" 1 


4.4X10" 1 


5.8X10- 1 


18) 


GT512mM10Bl 


(sccd3) 


4.0X10" 1 


4.5X10" 1 


4.6X10- 1 


5.3X10- 1 


3.9X10- 1 


4.2X10" 1 


5.3X10- 1 


19) 


GT256mM10B3 


(sccd3) 


3.4X10" 1 


5.1X10" 1 


1.6X10" 1 


2.6X10" 1 


1.8X10" 1 


2.0X10" 1 


3.1X10" 1 


20) 


GT512mM10B3 


(sccd3) 


2.9X10- 1 


5.0X10- 1 


1.4X10" 1 


2.5X10" 1 


1.7X10- 1 


1.9XKT 1 


3.0X10- 1 


21) 


GT256mM10B10 


(seeds) 


2.4X10" 1 


2.4x10+° 


2.8X10- 1 


1.6X10" 1 


3.5X10" 1 


3.6X10" 1 


2.3X10" 1 


22) 


GT128cM10 




1.6x10+° 


n/a 


n/a 


n/a 


n/a 


n/a 


n/a 


2.3) 


GT256cM10 




2.3x10+° 


2.5x10+° 


2.2x10+° 


1.1x10+° 


2.2x10+° 


2.3x10+° 


1.1x10+° 


24) 


GT512cM10 




2.8x10+° 


2.5x10+° 


2.2x10+° 


1.1x10+° 


2.2x10+° 


2.3x10+° 


1.1x10+° 


25) 


GT256sM20 




3.3X10- 1 


1.4X10- 1 


3.7X10- 1 


8.6X10- 1 


3.7X10- 1 


3.5X10- 1 


8.5X10- 1 


26) 


GT256mM20 




5.9X10" 1 


4.5X10" 1 


1.1x10+° 


1.0x10+° 


9.4X10" 1 


9.2X10" 1 


9.9X10- 1 


27) 


GT256cM20 




4.8x10+° 


2.3x10+° 


3.4x10+° 


2.0x10+° 


4.0x10+° 


4.1x10+° 


1.9x10+° 


28) 


GT256sM50 




3.8X10- 1 


l.lxlO -1 


9.3X10- 1 


1.8x10+° 


8.6X10- 1 


8.3X10" 1 


1.7x10+° 


29) 


GT512sM50 




4.4X10" 1 


9.9xl0~ 2 


8.8xl0 _1 


1.8x10+° 


8.2X10" 1 


7.9X10" 1 


1.7x10+° 


30) 


GT256mM50 




5.5X10- 1 


2.4X10" 1 


1.8x10+° 


2.0x10+° 


1.6x10+° 


1.5x10+° 


1.9x10+° 


31) 


GT512mM50 




6.8X10- 1 


2.3X10- 1 


1.7x10+° 


2.0x10+° 


1.5x10+° 


1.5x10+° 


1.9x10+° 


32) 


GT256cM50 




4.7x10+° 


1.9x10+° 


6.0x10+° 


3.9x10+° 


7.6x10+° 


7.7x10+° 


3.7x10+° 


33) 


GT512cM50 




7.3x10+° 


1.8x10+° 


6.1x10+° 


4.0x10+° 


7.7x10+° 


7.8x10+° 


3.7x10+° 


34) 


GT1024cM50 




9.1x10+° 


1.8x10+° 


6.1x10+° 


4.0x10+° 


7.8x10+° 


7.9x10+° 


3.8x10+° 



Notes. Column (1): simulation model. Column (2): SFRff measured in the simulations. Columns (3—8): Theoretical SFRff 
computed for the simulation parameters ct v i r , M, b, and /3 (or equivalently Ma) listed in Table J2] in the KM (3), PN (4), and HC 
(5) theories, as well as, in the multi-freefall KM (6), multi-freefall PN (7), and multi-freefall HC (8) theories, using the best-fit MHD 
parameters from Tablej3] No theoretical values were computed for the GT128sM10 and GT12 8cM10 simulations because they only used a 
numerical resolution ofT28 3 cells (see the discussion on numerical convergence in Section| 4.2| l. 

greater than zero, Egp > 0, and extract the correspond- 
ing pixel in the gas column density map, which gives 
E gas in units of M Q pc~ 2 for that pixel. The SFR col- 
umn density is computed by taking the sink particle col- 
umn density £sf of the same pixel and dividing it by 
a characteristic timescale for star formation, igp, which 
yields Ssfr = ^sfAsf in units of M Q yr _1 kpc~ 2 . The 
simplest choice for t^F is a fixed star formation time, 
tsF = 2 Myr, based on an estimate of the elapsed time 
betw een star formation and the end of the Class II phase 
(e-g 



the t oy 



Evans et al.]|2009[ [Covey et al.l|2010[ ). This is also 
f ado pte'crDy Lada et al. ( |2010J ) and Heiderman 



et al.| ( |2010[ ) to convert young stellar object (YSO) counts 
into an SFR column density, so we use it here as the stan- 
dard approach. However, we also experimented with two 
other choices for tsF and present a comparison of those 
choices below, all yielding similar results. 

The result of the proced ure explained above is plot- 
ted in panel (a) of Figure 11 It shows a scatter plot 
of Esfr versus S gas measurea in all the maps produced 
from our simulations listed in Table [2] (except for the two 



low-resolution, 128 3 -simulations) for a star formation ef- 
ficiency SFE = 1% (blue) and SFE = 10% (red). Thus, 
each pixel shown in panel (a) of Figure 11 is one pair of 
(S gas , Ssfr) extracted for each simulation and each pro- 
jection direction. By combining all data of maps from 
the three principal projections in x, y, and z, we in- 
crease the statistical sample for each model by about a 
factor of three on average. A total number of 3.5 x 10 3 
and 1.2 x 10 4 simulation pixels for SFE = 1% and 
SFE = 10%, respectively, contribute to the scatter plots 
in Figure fTTJ We also add contours of the (£ gas ,I]sFR) 
distribution, with two contour levels for SFE = 1% (blue 
contours) and SFE = 10% (red contours). The thick 
contours enclose 50% of all simulation pixels and the 
thin contours enclose 99%. The contours help to eas- 
ily identify the underlying probability distribution of the 
scattered data points. 

The simulation data have a broad probability distribu- 
tion with a clear positive correlation between SgpR and 
S gas . The data for SFE = 10% are shifted to higher Esfr 
and lower £ gas compared to the SFE = 1% distribution 
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Fig. 11. — (a): Star formation rate column density Ssfr vs. gas column density S gas measured in the GRAVTURB simulations listed in 
Table[2]for a star formation efficiency SFE = 1% (blue) and SFE = 10% (red), respectively. Two contour lines for each SFE are drawn. The 
thick contours enclose 50% of all (£ ga s, £sfr) simulation pairs, centered on the peak of the distribution, while the thin cont ours enclose 
99%. (b): Same as (a), but only the contours of the simulations are drawn, and observational data of Galactic clouds from [Heiderman] 
|et al | (2010] are superimposed. The individual data points are labeled in the legend of the bottom panels (Taurus: filled black box, (Jlass 1 
Y BUs and Flat YSOs: green and red stars and upper-limits shown as downward-pointing triangles, HCN(1— 0) Clumps: golden diamonds, 
and C2D+GB Clouds: dark blue boxes). The simulation data in panels (a) and (b) are plotted for a local core-formation efficiency e = 1, 
the value e xpe cted without any local feedback from YSOs. (c): Same as (b), but the simulation data were transformed to e = 0.5 using 
Equations }47| , which changes the GRAVTURB contour s co mpared to (a) and (b). The value e = 0.5 was determined by fitting the 
simulation data to the observational data using Equation \48\ , suggesting local efficiencies of e rj 0.3-0.7 for an assumed SFE « 1%-10% 
in the observed clouds, (d): Same as (c), but for the simulation maps smoothed to 4x coarser resolution, demonstrating the effect of 
observing the simulated clouds with reduced telescope resolution. 

because more gas is accreted by sink particles and thus 
removed from the gas phase at higher SFE. If we were 
to fit power laws to the distributions, the slopes would 
be in the range 1-2, i.e., Ssfr °c E ga s with somewhat 
flatter slopes at higher SFE. 

6.2. Galactic Observations o/£ gas and Ssfr 
To compare the simulation data with observations, we 



add data of Galactic clouds from Heiderman et al. (20101 
in panel (b) of Figure 11 superimposed on the simula- 
tion contours. The observational data are from Galactic 



Cores-to-Disks (C2D) and GB surveys dEvans et al. 2009) 


of massive dense clumps |Wu et al. 12010 


) , and ot the 


Taurus molecular cloud (.Pineda et ai.|20IU 


Rebull et al. 



2010 1. 

tours of panel (a) fall in the range of the observational 
data, however, the simulation data show slightly higher 
Ssfr than the observational data, on average. This is 
not surprising, given that our simulations did not include 
any local feedback from YSOs. It is known, however, 
that young stars eject a significant amount of accreted 
material, thereby reducing the overall accreti on rate due 
to feedback from jets, winds, and outflows (Wardle & 



Koenigll [19931 IKonigl & Pudritz| [20001 l Beuther et al - 
20021 IPudritz et al |2007||Peters et al.|2011||Seifried et al. 
2011ap . Hence, only a fraction e < 1 of the in-falling gas 
actually ends up on the protostar. 

The local core-formation efficiency is parameterized by 
the factor e in Equation frf\, from which all the SFRff- 
models in Section [2] were derived. Since there is no feed- 
back in our simulations, e — 1 by definition. However, 
we can devise a correction to account for e < 1. For this, 
we simply have to multiply the original Ssfr f° r e = 1 
by a given e < 1. To conserve mass, we also have to ac- 
count for the fact that a fraction (1 — e) was not accreted 
and remained in the gas phase due to local feedback. 
This means we have to increase E gas according to the 
reduction of Ssfr, such that S to t = S gas + Ssf with 
Esf = Ssfr^sf is conserved. Given our simulation data 
Sgas an d Ssfr with e = 1, we can compute values Sg FR 
and S' for e < 1 according to the following equations: 



s 'sfr( £ ) = £ S S fr , 
E' (e) = S 6a8 + (l- e ) S SF 



(47) 



Using these expressions, we can correct our simulation 
data to follow more realistic values of the local efficiency 
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Fig. 12. — Same as panel (c) i n F igure |11| but here we compute the simulation EgpR with two other methods, both different from the 
standard method used in Figure [ll] where igFR = Egp/(2 Myr). Left: EgpR = SsFAft(Po)i i- e -, the sink particle column density Egp is 
divided by the global freefall time at the mean density po of the simulation in which the Egp pixel was found. Right: EgpR = Egp/tff (E gas ), 
i.e., we divide Egp by the local freefall time of the gas for each pixel, %(Eg aa ) = ^/37rL/(32GE gas ) with the line of sight L of the 
corresponding simulation model. Both po and L are listed in Table ||] Some minor differences compared to panel (c) in Figure 11 
apparent, but the overall agreement between simulations and Galactic cloud observations remains good, irrespective of the method li 
to define EgpR in the simulations. 



arc 
used 



(see al so the discussion o f e in S ection 2.3) 



The Heiderman et al. 



SFR column 



(|2010j sampleof 
densities tor Galactic clouds shown in panel (b) of Fig- 
ure [TT] is rather broad and presumably covers different 
evolutionary stages of the clouds, such that a single SFE 
for the whole sample is quite unlikely. However, since we 
are currently lacking additional information about the 
SFE in the observational sample, we can reasonably as- 
sume SFEs in the range 1%— 10% in the observ ational 
data ( |Evans et al |2009| |Federrath & Klessen|2012[ Paper 
II). In order to find the best-tit local efficiency parame- 
ter e, we fit our simulated distribution p s i m (E gas , Eg FR ) 
to the observe d di stribution Pobs(E gas , Esfr), by apply- 
ing Equations (47 1. To do this, we compute the sum of 



the squared differences A 2 between the two distributions, 
which have both been sampled to the same (E gas , Esfr) 
grid with indexes i, 



A 



[psim(Eg as ,i, E' SFRi i) — p bs(E gaSi i, Esfr,;)] , 

(48) 

for SFE = 1%, 3%, and 10% and for 21 local efficiencies, 
e = [0, 1] in steps of de = 0.05. For each given SFE, 
we search for the minimum of A 2 as a function of e. 
This procedure yields best-fit values of the local efficiency 
parameter e = 0.7, 0.5, and 0.3 for SFE = 1%, 3%, and 
10%, res pectively, in our compar ison of simulation data 
with the Heiderman et al. (2010) Galactic cloud sample. 
The simulation data modified to a local efficiency of 



Lada et al.| ( |2010[ ) and |Gutermuth et al.| ( |2011[ ), showing 
that Esfr can vary by more than an order of magnitude 
at any given E gas . 

Considering the uncertainties in the SFE from the ob- 
servations and the uncertainties in the simulations, the 
overall agreement is encouraging. The HCN(l-O) obser- 
vational data points of molecular clumps are at the lower 
end of the distribution, but are still consistent with the 
simulation data. Possibly, the molecular clumps have 
a systematically smaller SFE because they are larger 
structures compared to the YSOs, such that the molecu- 
lar clumps fall slightly below the general trend. How- 
ever, this can only be tested when estima tes of the 
cloud SFEs beco me available (see Paper II, Federrath 
& Klessen 2012 1. The Taurus data point as well as a 
few ot the YSU data in the range log 10 E gas s» 1.4-2.8 
also lie at the Iow-Esfr end of the distributions ob- 
tained in the simulations. This might be caused by an 
enhanced magnetic-field influence for these objects. For 
instance, Taurus seems to be trans-Alfvenic rather than 



11 



e = 0.5 are sho wn in panel (c) o f Figu re 
with the original [Heiderman et al. (2010) data. Assum 
ing that the observational data have an SFE between 1% 
and 10%, the local efficiency parameters would be in be- 
tween e = 0.3 and 0.7. This is in good agreemen t with 
theoretical models for e ( Matzner k McKee||2000[ ) , with 
numerical simulations incl uding outflow feedback ( jWang 
et al.p OlO; Scifried et al. 20121), and with observational 
estimates (Bcu thcr et al.p002[ and the discussion on e 
in Section 2.3) 



We note that the simulation data in Figure 11 are fur 



thermore consistent with the Galactic cloud samples in 



super- Alfvenic ( |Heyer k Brunt 2012 1, leading to a re- 
duced Esfr as discussed in Section |4.1| Only one of our 
MHD simulations approaches this strongly magnetized 
regime (GT256mM10B10 with M A ~ 1.3; see Table 
where anisotropies induced by the magnetic field start to 
become important. 

6.3. Influence of Telescope Resolution and Choice o/^sf 

We test the effects of telescope beam smoothing in 
panel (d) of Figure |TTj Panel (d) is identical to panel (c), 
except that the simulation data were smoothed to grids 
with resolution iV rcs /32, i.e., four times coarser resolu- 
tion compared to the contours shown in panel (c). The 
increased beam smoothing results in distributions with 
somewhat smaller E gas and Esfr, best seen by compar- 
ing the positions of the thickest contours between panels 
(c) and (d). However, the overall agreement of the simu- 
lation data with the Galactic cloud sample is still good, 
even when the resolution is decreased by a factor of four. 

In Figure 12 , we study the influence of different choices 
for the star formation timescale t^F- The two panels 
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are identical to panel (c) in Figure 11 except for the 
method by which Esfr = EsfAsf was com P u ted in the 
simulations. The left panel adopts isF = ^ff(po), i.e., 
the sink particle column density Esf is divided by the 
freefall time at the mean density p of the simulation in 
which the Esf pixel was found. In the right panel, we 
use tsF = iff(£ gas ) = \Z3nL/ (32GE gas ), i.e., instead of 
taking the global mean free-fall time, we take the local 
freefall time of the gas in each pixel. The contours differ 
slightly between those two last choices and between our 
standard choice of fixed isF = 2 Myr in Figure[lT] but the 
overall agreement between simulation data and Galactic 
observations is similar in all three cases. 

6.4. Comparison with Extragalactic Measurements 

Figures [TT] and [12] indicate some power-law correla- 
tion of the form Esfr oc Eg^ s (albeit with significant 
scatter) , similar in exponents N » 1-2 to the Kennicutt- 



Schmidt relation (Schmidt 1959 Kennicutt 1998) and 
follow-up measurements for molecular gas (e.g., |Won; 



fe Blitz|[2002l |Gao fe Solomon| [20041 [Bigiel et al.|| 2HB8 



Kennicutt & Evans 2012[ ). However, the measured val- 
ues of Usfr m our numerical sample are larger than the 
extragalactic values of Esfr and larger than theoretical 



estimates for that regime (e.g., Krumholz et al 
by about 1-2 orders o f magnitude 



2009 1 

i'he Galactic mea- 



surements of Esfr by Heiderman et al. (2010) in Fig- 



11 



and by Lada et 



values of E S fr that are 



t al j |2010 

e 1-2 ordei 



however, also show 
rs of magnitude above 



the extragalactic measure ments with a scatter of ab out 
1-2 orders of magnitude. Heiderman et al. (2010) ex- 
plain this difference between Galactic and extragalactic 
measurements of Esfr with the different telescope reso- 
lutions available for both regimes and thus the different 
areas over which the measurements of E gas and Esfr are 
averaged. Both disk-averaged and spatially-resolved ex- 
tragalactic measurements only provide highly smoothed 
images, mixing both star-forming and non-star-forming 
gas. Taking these factor s into a ccount and correcting for 
them, Heiderman et al.] ( |2010[ ) conclude that the extra- 
galactic (iJg as , Usfr] relations are in agreement with the 
Galactic measurements. Indeed, decreased telescope res- 
olution (or equivalently observing a region at greater dis- 
tance) reduces Esfr, but also E gas , as demons trated here 
by comparin g panels (c) and (d) of Figure 11 |KrumhoTz] 
et al.| (2012) argue that both Galactic ana extragalac- 
tic measurements are consistent with a local star forma- 
tion law, correlating Esfr with E gas /isF, where isF "is 
the freefall time evaluated at the density averaged over 
length scales comparable to the outer scale of turbulence, 
regardless of the mean density of the region being ob- 
served" . This seems to be a rather especial definition. 
Our experiments with three different definitions of igF i n 
Figures [TT] and [12] do not exclude or prefer any particular 
choice for tsF inuie Galactic cloud sample studied here. 
After acceptance of this work, we also learned about a 
recently submitted paper on a theoretical mod el for the 
EsFR-to-Eg as relation by Renaud et al. ( 2012[ ), which is 
consistent with our findings tor Galactic clouds, favoring 
a non-universal behavior of the star formation relation. 

The sim ulations and the observational data shown in 
Figure 11 are generally in very good agreement. The 
variations of the observed SFRs in different clouds by 



up to two orders of magnitude for a given value of E gas 
(jMooney fc Solomon | |1988| |Lada et aL]|2010| [Heiderman 



et al. 2010p and the different sca ling relations of Esfb. 



versus E gas (Suzuki et al. 2010) might thus be a re 



suit of different physical conditions in Galactic as well 
as extragalactic molecular clouds. As shown above, star 
formation is primarily controlled by the forcing and the 
sonic Mach number of the turbulence, with the magnetic 
field having a secondary effect. Molecular clouds cover a 
range of values for these physical parameters and differ- 
ent combinations of those, providing an explanation for 
the observed scatter in SFRs. 

7. DISCUSSION AND LIMITATIONS 

Here we discuss limitations of the analytic theories for 
the SFR from Section [2] the numerical simulations from 
Sections [3H51 and limitations of the comparison of both 
theory and simulations with observations in Section [6] 



7.1.1. 



7.1. Analytic Theories 
Non-log-normal Effects in the Density PDF 



One limitation of the current analytic theories for 
SFRff is the assumption of a perfect log-normal PDF 
of the gas density, Equation (nl), in the derivation of 
the SFR integral, Equation (CT, which affects all six 
analytic theories (Table [lj similarly. Even though a 
log-norma l PDF is expected for p urely isothermal tur- 
bulence ( |Vazquez-Semadeni 1994), intermitten cy intro- 
duces skewness and kurtosis in the distributions ( Klessen 



2000 



Kritsuk et al. 2007] IBurkhart et al. 2009 , which 



becomes st ronger for more compres sive forcing ([Schmidt 
et al.||2009| iFederrath et al.|2010b| and for higher Mach 
numbers ((Konstan din et al.||2012b| . Temperature varia- 
tions can also introduce deviations from perfect log nor- 
mals in the wings of the distributions. This occurs, for 
instance, if the turbulence is modeled with a polytropic 



unity (|Passot & Vazquez-Semadeni 


1998 


Jappsen et al. 


20U5). However, wJ 


len a 



Li et al. 
detailed. 



ing model is used instead of a polytropic EOS, the PDF 
of the main molecular gas co mponent, H 2 , follows a log- 
normal distributio n very well ([Glover fc Mac Low|2007a[ 
Glover et al.||2010[|Shetty et al.||201l[ |Micic et al.J|201; 



I'he strongest deviations from log-normal PDF arise 
when the gas starts to collapse due to se lf-gravity, pro 
during power-law tai l s at high densitie s 



Dib fc Burkert][2005l |FcdciTathlt~aT]|20(]8a[ | Vazquez 



2000 



Semadeni et ai.||20081 |Gho fc Kim||20il[ |Kritsuk et al 
2011al|Ballesteros-Paredes et al.|20111|Collinset al.|2012 



Safranek-Shrader et al. 2012[ ), which has been observec 
in the column density PDFs of clouds that have already 



forme d stars (Kainulainen et al. 



2012|. One might thu s argue that star formatioh~ might 



2009 Schneider et al. 



accelerate over time ( Cho & Kim |2011[ Collins et al. 



2012|. In our numerical experiments, we see that after 



an initial transient acceleration of SFE(£) in Figures [7] [9] 
the SFR becomes fairly constant in most of the numeri- 
cal models for SFE > 4%. This taken together with the 
good fit-quality of SFRff (theory) to SFRff (simulation) 
obtained for the multi-freefall KM and PN models in Fig- 
ure [10] suggests that the development of power-law tails 
in the density PDF during star formation does not sig- 
nificantly affect star formation itself. Using a log-normal 
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PDF in the analytic theories to estimate SFRg seems 
to be a reasonably good approximation. From a cer- 
tain perspective, we could say that the initial conditions 
for star formation are basically determined by the log- 
normal part of the PDF. In regions that form stars, the 
PDF develops a power-law tail at high densities, which 
is a result (or a byproduct) of star formation, but does 
not necessarily affect the process of stellar birth itself. 
We d iscuss this further in Paper II ( Federrath & Klessen 
2012 1, where we present the density FDFs of the Simula- 



tions, showing the development of power-law tails when 
star formation sets in, consistent with the assumption 
that the power-law tails obse rved in molecular clouds 



correlate with star fo rmation (Kainulainen et al. 
Schneider et~aL|2012| ). 



2009 



7.1.2. Anisotropies in Sub-Alfvenic Turbulence 

The present analytic theories only work for super 
Alfvenic turbulence because Equation ^ and ^ break 
down for Ma ;$ 2 (see the discussion m Molina et al. 
2012|. All theories assume statistical isotropy, which is 
only fulfilled in the trans- to super-Alfvenic regime of 
turbulence studied here. 

7.1.3. Vmal Parameter 

The virial parameter in Equation ( |15[ ) only applies 
to spherical, uniform-density clouds. Tn the compari- 
son with numerical simulations (columns 10 and 11 in 
Table p2|), it became clear th at t he virial parameter, 



^kin/ 



Egr&v], Equation (16), based on the spatial 



gas distribution can be more than an order of magni- 
tude differ ent from the virial parameter estimated by 
Equation (15). This is because turbulent interstellar 



gas is concentrated in fractal-like structures that differ 
significantly between solenoidal and compressive forc- 
ings, and between different sonic Mach numbers (see 
Figure [6]), even when the total mass is identical. How- 
ever, we also tested using a v ir,o instead of a v i r in the 
theory-simulation comparison of Section [3J Doing so 
yielded similar fits to the ones shown in Figure [10] and 
listed in Table [3j yet with somewhat larger x^ed m some 
cases. We thus preferred to use the direct computation 
of a v ir in the simulations, Equation (16), which provides 
a more meaningful description of the dynamical state of 
the clouds. In the derivation of the analytic models in 
Section [2j how ever, we use the simple definition given by 
Equation (fl5| because it can be treated analytically. 



7.2. MHD Simulations 

7.2.1. Approximation of SFRg as Constant Over Time 

In both the theory and MHD simulations, we approx- 
imate SFRg as constant over time. Figures [7|j9] show 
that this is a reasonable assumption for SFE > 4%, but 
the initial acceleration of SFRg when SFE < 4% is more 
complicated and is not accounted for in the present the- 
ory and simulations. In real molecular clouds, the SFRg 
might also change over time, depending on the evolution- 
ary stage of a cloud, or on environmental parameters. 

7.2.2. Limited Numerical Resolution 

Our numerical resolution studies in Figures [7] and [10] 
show that SFRg converges with increasing resolution in 
the numerical simulations. However, some models still 



differ by a factor of 2-3 from the best analytic predic- 
tions. In particular, the very high Mach number simula- 
tions with A4 ~ 50 are not converged at a resolution of 
256 3 and only marginally resolved with 512 3 cells. How- 
ever, the 1024 3 -simulation GT1024cM50 with compres- 
sive forcing at M. ~ 50 seems reasonably well converged 
as suggested by Figure [l0| (model #34). The lower-Mach 
number simulations typically agree within a factor of 
1.5 with the best analytic theories (see Table E|, which 
is similar to the typical statistical variation induced by 
different random realizations of the turbulence (see the 
comparison of three different random seeds in Figure [7]). 

7.2.3. Periodic Boundary Conditions 

Our numerical simulations are highly idealized in that 
the boundary conditions are periodic. Real molecu- 
lar clouds are embedded in the larger-scale interstel- 
lar medium and eventually in galaxies, which sets their 
boundary conditions. Our choice of boundaries intro- 
duces some uncertainties, e.g., in the virial parameter 
because the gravitational energy E'grav entering a v - lr de- 
pends on the choice of boundary condition. The other 
extreme would be to in itialize a cloud in isolation as done 
in related studies (e.g., Bate et al.|200 3||Clark et al.|2005 
Krumholz et al.||2007| |Price & Bate 



2008, 2009; Smit 



et al.| 2008; Federrath et al. 2010a Walch et al. 2010 



15; 
th 



Girichidis et al.|201ip . This is similarly artificial because 



real clouds are not isolated, but exist in a large-scale 
interstellar web of filaments and other clouds. 

Here, we test the analytic theories introduced in Sec- 
tion [2] with s uch simulations o f isolat ed star formation. 
For instance, Girichidis et al. (2011) modeled isolated 
clouds with different density profiles and an initial turbu- 
lent perturbation, i.e., impulsive turbulent forcing. Since 
the clouds with initial power-law or Bonnor-Ebert pro- 
files already assume a stage of previous evolution that 
may have led to such a density profile, we prefer to com- 
pare the more basic, simple init ial condition when the 
density field is initially uniform. Girichidis et al. (20111 



modeled such a uniform density profile with a mixed {b 
0.4) turbulent perturbation with two different random 
seeds, in which the sonic Mach number M. — 3.3 for their 
simulation TH-m-1 and M = 3.6 for TH-m-2. The sim- 
ulations did not include magnetic fields, so /3 — > oo . The 



virial param eters are in the range a v ir = 1 _ 2 ( Girichidis 
et al.||2012|), depending on the time interval and spa- 



tial range chosen to determine a v - lT , which exhibits some 
temporal and spatial variation. Using the best-fit multi- 
freefall PN parameters determined from Figure [10] and 
Table [3] (l/<f> t = 0.47 ± 0.16 and 6 = 1.0 ± 03), an 
average virial parameter a. y - a — 1.5, an average Mach 
number of A4 = 3.45, and b — 0.4 for mixed turbu- 
lence, we find SFRg(multi-ff PN) = 0.56 b y e valuating 
Equation (41) with s cr jt.PN from Equation (30). Taking 
the uncertainties in the fit parameters l/(f>t and 9, as 
well as the uncertainty in a V i r = 1 _ 2 and M. = 3.3-3.6 
into account, we find the analytic multi-ff PN predic- 
tion SFRg(multi-ff PN) = 0. 56 ± 0.35 for both TH- m-1 
and TH-m-2 simulations by Girichidis et al. (20111. A 



very similar prediction is obtained using the multi-treefall 
KM model instead of the multi-freefall PN model with 
the corresponding parameters listed in Table [3j From 
a linear fit to the evolution of the total accreted mass 
versus time in the TH-m-1 and TH-m-2 simulations, we 



30 



Fcdcrrath & Klessen 



find SFR ff (TH-m-l) » 0.67 and SFR ff (TH-m-2) « 0.61 
(Girichidis et al. 2011 Figure 7), in very good agree- 
ment with the analytic model prediction, indicating that 
different boundary conditions do not severely affect our 
results and conclusions concerning SFRff. 

7.3. Observations 
Assuming a uniform SF E = 1%-10% in the obse rved 



Galactic cloud sample by Heiderman et al. (20101, we 
estimated the local core- format ion efficiency parameter 
e = 0.3-0.7 with the best-fit value e ~ 0.5, by fitting 
our numerical simulations to the observed distribution 
in Figure [TT] There are three major uncertainties in this 
comparison of the simulations and observations. 

First, the SFE in the observed sample is not known. 
We reasonably assumed SFE = 1%-10%, but some of the 
individual clouds may not fall in this range. Moreover, 
there could be a systematic correlation of SFE with gas 
column density S gas , which is not accounted for. For 
instance, the HCN(l-O) molecular clump data shown in 
Figure [TT] has potentially smaller SFE on average than 
the YS(J da ta because smaller scales tend to exhibit 
higher SFE ( |McKee fc Ostriker| |2007[ ) . For instance, 
it seems plausible that SFE approaches the local core- 
efficiency e, once scales as small as a single core are con- 
sidered. In contrast, giant molecular cloud complexes as 
a whole typica lly only have SFEs of a few percent at most 
(see Paper II, Federrath & Klessen 2012). 

The second uncertainty is the effect of the telescope 
resolution. Lower resolution (or observation of a very 
distant region, e.g., a whole galaxy) inevitably means 
that the observed star-forming regions are smoothed over 
larger areas compared to a high-resolution observation of 
the same region. The effect of reducing the beam reso- 
lution by a factor of four in our synthetic observations 
of the simulated clouds is dem onstrated by comparing 
panels (c) and (d) in Figure 11 resulting in a relatively 



weak, but noticeable reduction of Ssfr and S gas . 

The third major uncertainty is the star formation 
timescale isF used to convert a given star formation col- 
um n de nsity Esf into a rate Esfr = ^sf / tsF ■ In Fig- 
ure [TTJ we adopte d a fixed tsF = 2Myr as often used 
by oE servers (e.g., Heiderman et al. |2010| Lada et al. 
20101. However, we studied two additional choices of tsF 
in Figure 12 one where tsF = ts{po) (division by the 
global frceiali time) and the other where isF = te(S gas ) 
(division by the local freefall time). Comparing these 
three choices for igp , we find that the resulting SgpR-to- 
S gas correlations change slightly, but the overall effect 
is rather weak. Given the broad distributions in both 



the simulation data and in the Heiderman et al. (20101 
Galactic cloud sample, it is hard to decide which method 
provides better agreement. They all seem to agree rea- 
sonably well within the observational range of Galactic 
clouds. 

Finally, we note a fundamental difficulty of estimat- 
ing actual SFRs or SFRg in observations. Cloud obser- 
vations are inevitably limited to a nearly instantaneous 
snapshot of the state of a cloud with respect to the rele- 
vant timescales for star formation, which exceed the life- 
time of a human being by orders of magnitude. How- 
ever, measuring a real SFR requires knowledge about 
the time evolution of the cloud, which is thus not avail- 
able. Strictly speaking, a direct measurement of the time 



derivative of star formation, i.e., the SFR is thus impos- 
sible in observations. This is why we can only make 
meaningful comparisons of star formation in simulations 
and observations based on th e m etho ds explained and 
applied in Section [6] (Figures nuand[l2|), but not the 
actual SFRs computed from the t 
formation. 



time evolution of star 



8. SUMMARY AND CONCLUSIONS 

We investigated the role of turbulence and magnetic 
fields for the SFR in molecular clouds. We compared 
theoretical models for the SFR with a comprehensive set 
of numerical magnetohydrodynamic simulations of core 
and star formation, and with observations of Galactic 
clouds. The main conclusion from this study is that the 
SFR depends on four parameters: (1) the virial parame- 
ter, a v ir = 2£'kin/|-Egrav|; (2) the sonic Mach number M; 
(3) the turbulent forcing parameter b (solenoidal, mixed, 
compressive); and (4) the strength of magnetic fields, pa- 
rameterized by plasma (3 = 2A4 A /J\A 2 with the Alfven 
Mach number Ma- 

Our simulations are in good agreement with SFR 
column densities and gas column densities of observed 
molecular clouds. We suggest that variations of the four 
basic, dimensionless parameters can explain the scatter 
in the observations. Given that molecular clouds seem 
to have an a V i r of order unity, the most important pa- 
rameters controlling the SFR are the sonic Mach number 
M. and the turbulent forcing of a molecular cloud, with 
magnetic field having a secondary effect. The turbulent 
forcing can be parameterized by b in Equation Q . It is a 
measure for the fraction of energy excited in the form of 
compressive modes in a turbulent cloud. We distinguish 
solenoidal (divergence-free) forcing (6 = 1/3) from com- 
pressive (curl-free) forcing (b = 1), as well as mixtures of 
both (1/3 < b < 1). We find that the SFR decreases with 
increasing magnetic pressure, but only by a factor of two. 
The sonic Mach number can change the SFR by a factor 
of 4-5, while b can introduce order-of-magnitudc differ- 
ences in the SFR, emphasizing the role of the turbulent 
forcing for star formation. 

8.1. Analytic Theories for SFRff 

1. In Section [2j we derived six analytic models for 
the SFR per freefa l l time , SF Rg: the original 



Kmmholz & McKee 



(120051 KM),|Padoan fc Norcfl 
lund| ( |20lTj PJN), anJ lHennebelle fc Chabrier|(|2mTr 

HU) models and the multi-freefall KM, Pi\, anc 
HC models, which are all based on an integral over 
the density PDF, Equation (fTl), leading to differ- 
ent analytic solutions for SFRff, summarized in 
Table [TJ They all yield a dimensionless SFR per 
freefalltime, SFRff, based on Equation Q, which 
can be transformed to a real SFR with units of 
M Q yr _1 by applying Equation M^. 

2. We extended the (multi-freefall) KM and (multi- 
freefall) HC theories to include magnetic fields by 
introducing a ma gnetic-pressure correction given 



by Equation (17), which allows us to replace the 
sound speed by an effective magnetosonic speed 
given by Equation ( 18 ) or ( 19 1 for super- Alfvenic, 



isothermal turbulence 
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3. We analyzed the basic dependencies of all six the- 
ories on the four parameters listed above. SFRff 
decreases with increasing virial parameter a v i n 
while it increases with increasing sonic Mach num- 
ber A4 in the best multi-freefall theories (see Fig- 
ure [I]). Varying the forcing parameter b from 
purely solenoidal forcing (6 = 1/3) to purely com- 
pressive forcing (b = 1) leads to a higher SFRff 
by more than an order of magnitude (Figure [2]) . 
Stronger magnetic fields parameterized by decreas- 
ing plasma /3 (or equivalently decreasing Alfven 
Mach number A^a) lead to decreasing SFRff (Fig- 
ure [3]). 

8.2. Numerical Simulations 

1. In Sect ions [3] and |4j we performed a set of numerical 
experiments of star formation, covering molecular 
cloud sizes and masses in the range L — 0.3 to 
200 pc and M c = 300 to 4 x 10 6 M Q (see Table [2) 
with solenoidal, mixed, and com pressive forcings 
of the turbulence (see Section 3.2 for details of the 
forcing) to test the analytic models. We also ran 
super- Alfvenic simulations with varying magnetic- 
field strength to test the influence of magnetic fields 
on the SFR. All simulations include sink particles 
to model core and star formation, allowing us to 
measure SFRff, depending on a v - n -, M, b, and /3. 

2. We computed the virial parameter a V i r = 
2-Eidn/l-E'gravl based on the uniform-density, spher- 
ical approximation given by Equation (15), and 
based on the actual, three-dimensional, mhomo- 
geneous gas distribution in the simulations. De- 
pending on the forcing and Mach number of the 
turbulence, we find that these two definitions can 
differ by more than an order of magnitude (com- 
pare columns 10 and 11 in Table [2|, which means 
that theoretical and observational estimates of a v j r 
based on a uniform-density, spherical approxima- 
tion must be considered with caution. 

3. The SFR converges wi th i ncreasing numerical res- 
olution (Figures u\ and 10). The statistical uncer- 
tainty in SFRff is about a factor of 1.5, indicated 
by comparing three different random realizations of 
the same parameter set (Figure [7]) , similar to the 
uncertainty introduced by limitea numerical reso- 
lution. 

4. We found that for our models with M. ~ 10, com- 
pressive forcing yields SFRs at least an order of 
magnitude higher than solenoidal forcing, empha- 
sizing the role of different turbulent energy injec- 
tion mechanisms for the SFR (Figure[7]). The cloud 
morphology also depends strongly on the type of 
forcing and sonic Mach number (see Figure pi). The 
SFR increases by a factor of about four for com- 
pressive forcing between A4 — 3 and A4 — 50 (Fig- 
ure [8]). 

5. Including magnetic fields in simulations with A4 ~ 
10 and mixed turbulent forcing, we found that the 
magnetic field is amplified in regions of core and 
cluster formation (Figure [5]), reducing the SFRff by 



a factor of two between purely hydrodynamic tur- 
bulence (Ma ^ 00) and trans- Alfvenic turbulence 
with A4a ~ 1-3 (see Figure[9]). This is a relatively 
small change in SFRff for such a fairly strong mag- 
netic field, compared to the dependence of SFRff on 
a v ; r , M, and b. However, magnetic fields do affect 
the morphology of the clouds even on large scales, 
and they reduce fragmentation (see Figure El), thus 
potentially having an important impact on the core 
and stellar IMF. 

6. A detailed comparison o f S FRff (simulation) with 
SFRff (theory) in Figure 10 showed that the multi- 



freefall analytic theories are generally better than 
the non- multi-freefall theories. The multi-ff KM 
and multi-ff PN models give the best fits to our 
simulation data (see Tables [3] and EJ) with reason- 
able best-fit model parameters, l/4>t ~ 0-5 for both 
multi-ff KM and PN models, as well as <f> x f=s 0.17 
for the multi-ff KM model, and 8 s» 1 for the multi- 
ff PN model, suggesting a close connection between 
the magnetothermal Jeans scale and the magne- 
tosonic scale, as well as turbulence driven on the 
outer, largest scales of molecular clouds. 

7. All numerical simulations agree with the multi-ff 
KM and PN theories within a factor of three, and 
come closer to the analytic prediction with increas- 
ing numerical resolution. This is an encouraging 
agreement, given that the modeled SFRs vary over 
two orders of m agni tude in our numerical simula- 
tions (see Figure 10). 



8.3. Comparison with Observations 

1. We compared our numerical simulations with ob- 
servations of the SFR column density Esfr as a 
function of the gas column density S gas , mea sured 



m Galactic clouds in Section [6[ (Figure 11). We 
showed that the simulations slightly overestimate 
the SFR compared to the observed clouds because 
we did not include any local radiative and mechani- 
cal feedback from young stellar objects, and hence, 
the local efficiency parameter e = 1 in our simu- 
lations, by definition. However, assuming a con- 
stant, global star formation efficiency in th e ob 
served clouds of SFE w 1%-10% (see Paper II, [Fed 
errath & Klessen|2012|), we can adjust our numeri 



cal simulation data with Equations ( 47 1 to account 
for e < 1. Doing so, we found the besWit local effi- 
ciency e « 0.5 (e = 0.7, 0.5, and 0.3 for SFE = 1%, 
3%, and 10%, respectively) for the observed Galac- 
tic clouds, which is in good agreement with theo- 
retical expectations, independent numerical simu- 
lations, and observations of individual protostellar 
cores. 

2. We studied the effects of telescope beam smooth- 
ing in panels (c) and (d) of Figure 11 and the ef- 
fect of varying the definition of the star fo rma tion 
timescale t$p to determine Ssfr in Figure [12] We 
found that both the telescope beam resolution and 
the definition of isF introduce minor uncertainties 
in our comparison between simulations and obser- 
vations. 
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3. The correlation between £„ 



ure 
Ssfr 



gas and E S fr in Fig- 
11 is consistent with power laws of the form 
oc E^, s with exponents N = 1-2 (albeit with 



, as with exponents N 
significant scatter), which is similar to extragalac- 
tic measurements of E gas -EsFR correlations. 

The overall agreement bet ween theory, simulations and 

is encouraging, 



observations in Figures 10 and 11 



sidering the simplifications inherent in the theoretical 
models, the limitations of the numerical simulations, and 
the uncertainties in the SFEs of the observed sample of 
clouds (see Section [7]). We conclude that supersonic, 
magnetized turbulence is a key process, likely control- 
ling the SFR of molecular clouds in the Milky Way and 
potentially in other galaxies. 
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